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Abstract 



A rigorous QCD analysis of the inclusive annihilation decay rates of heavy 
quarkonium states is presented. The effective-field-theory framework of non- 
relativistic QCD is used to separate the short-distance scale of annihilation, 
which is set by the heavy quark mass M, from the longer-distance scales asso- 
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dated with quarkonium structure. The annihilation decay rates are expressed 
in terms of nonpcrturbative matrix elements of 4-fermion operators in non- 
relativistic QCD, with coefficients that can be computed using perturbation 
theory in the coupling constant as{M). The matrix elements are organized 
into a hierarchy according to their scaling with v, the typical velocity of the 
heavy quark. An analogous factorization formalism is developed for the pro- 
duction cross sections of heavy quarkonium in processes involving momentum 
transfers of order M or larger. The factorization formulas are applied to the 
annihilation decay rates and production cross sections of S-wave states, up 
to corrections of relative order v^, and of P-wave states, up to corrections of 
relative order v'^. 
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I. INTRODUCTION 

Calculations of the decay rates of heavy quarkonium states into light hadrons and into 
photons and lepton pairs are among the earliest applications of perturbative quantum chro- 
modynamics (QCD) [|^,^^Q. In these early analyses, it was assumed that the decay rate 
of the meson factored into a short-distance part that is related to the annihilation rate of 
the heavy quark and antiquark, and a long-distance factor containing all the nonpertur- 
bative effects of QCD. The short-distance factor was calculated in terms of the running 
coupling constant as{M) of QCD, evaluated at the scale of the heavy-quark mass M, while 
the long-distance factor was expressed in terms of the meson's nonrelativistic wavefunction, 
or its derivatives, evaluated at the origin. In the case of S-waves P,^ and in the case of 
P-wave decays into photons , the factorization assumption was supported by explicit cal- 
culations at next-to- leading order in a^. However, no general argument was advanced for 
its validity in higher orders of perturbation theory. In the case of P-wave decays into light 
hadrons, the factorization is spoiled by logarithmic infrared divergences that appear in the 
QQ annihilation rates at order a^ [133 • Logarithmic infrared divergences also appear in 
relativistic corrections to the annihilation decays of S-wave states ||^. These divergences 
cast a shadow over applications of perturbative QCD to the calculation of annihilation rates 
of heavy quarkonium states. 

In this paper, we present a rigorous QCD analysis of the annihilation decays of heavy 
quarkonium. We derive a general factorization formula for the annihilation rates of S-wave, 
P-wave, and higher orbital-angular-momentum states, which includes not only perturbative 
corrections to all orders in a^, but relativistic corrections as well. Factorization occurs in 
the annihilation decay rates because the heavy quark and antiquark can annihilate only 
when they are within a distance of order 1/M, where M is the heavy-quark mass. Since, 
in the meson rest frame, the heavy quark and antiquark are nonrelativistic, with typical 
velocities f <^ 1, this distance is much smaller than the size of the meson, which is of order 
l/{Mv). Factorization involves separating the relativistic physics of annihilation (which 



involves momenta p ~ M) from the nonrelativistic physics of quarkonium structure (which 
involves p ~ Mv). A particularly elegant approach for separating relativistic from nonrela- 
tivistic scales is to recast the analysis in terms of nonrelativistic quantum chromodynamics 
(NRQCD) 0, an effective field theory designed precisely for this purpose. NRQCD con- 
sists of a nonrelativistic Schrodinger field theory for the heavy quark and antiquark that 
is coupled to the usual relativistic field theory for light quarks and gluons. The theory is 
made precisely equivalent to full QCD through the addition of local interactions that sys- 
tematically incorporate relativistic corrections through any given order in the heavy-quark 
velocity v. It is an effective field theory, with a finite ultraviolet cutoff of order M that 
excludes relativistic states — states that are poorly described by nonrelativistic dynamics. 
A heavy quark in the meson can fluctuate into a relativistic state, but these fluctuations are 
necessarily short-lived. This means that the effects of the excluded relativistic states can 
be mimicked by local interactions and can, therefore, be incorporated into NRQCD through 
renormalizations of its infinitely many coupling constants. Thus, nonrelativistic physics is 
correctly described by the nonperturbative dynamics of NRQCD, while all relativistic effects 
are absorbed into coupling constants that can be computed as perturbation series in as{M). 
The main advantage offered by NRQCD over ordinary QCD in this context is that it is 
easier to separate contributions of different orders in v in NRQCD. Thus, we are able not 
only to organize calculations to all orders in a^, but also to elaborate systematically the rel- 
ativistic corrections to the conventional formulas. Furthermore, we provide nonperturbative 
definitions of the long-distance factors in terms of matrix elements of NRQCD, making it 
possible to evaluate them in numerical lattice calculations. Analyzing S-wave decays within 
this framework, we recover, up to corrections of relative order v"^, the standard factorization 
formulas, which contain a single nonperturbative parameter. At relative order f ^, the decay 
rates satisfy a more general factorization formula, which contains two additional indepen- 
dent nonperturbative matrix elements. Our results for P-wave decays into light hadrons are 
even more striking, as we have discussed in Ref. [0. Up to corrections of relative order v"^, 
the factorization formula for these decay rates is the sum of two terms. In addition to the 



conventional term, which takes into account the annihilation of the QQ pair from a color- 
singlet P-wave state, there is a second term that involves annihilation from a color-octet 
S-wave state. The infrared divergences encountered in previous calculations are absorbed 
into the matrix element of the color-octet term. 

Our presentation is organized as follows. In Section 0, we first review NRQCD in gen- 
eral, emphasizing the velocity-scaling rules, which are used in separating contributions of 
different orders in v. We then discuss the space-time structure of the annihilation of heavy 
quarks and antiquarks and explain how the effects of annihilation can be taken into account 



in NRQCD by adding local 4-fermion operators to the effective lagrangian. In Section |T|, 
we analyze the matrix elements of the 4-fermion operators. We discuss their scaling with 
f , the constraints on them that follow from heavy-quark spin symmetry, their relations to 
Coulomb-gauge wavef unctions, and their dependences on the factorization scale. In Sec- 



tion 1^ we apply our formalism to the annihilation decays of S-wave quarkonium states, up 
to corrections of relative order f ^, and to P-wave decays, up to corrections of relative order 
f^. In Section 0, we sketch the derivation of our results in a more conventional perturbative 
approach to factorization. In Section [V|, we develop an analogous factorization formalism 
for calculating the production cross sections of heavy quarkonium. In the concluding section, 
we compare our formalism with previous approaches to the annihilation and production of 
heavy quarkonium, and we summarize the current status of calculations of annihilation and 
production rates. 



II. NRQCD 

We begin this section with a brief discussion of the various momentum scales involved 
in heavy quarkonia. Nonrelativistic QCD (NRQCD) |^ is our major tool for resolving the 
different momentum scales involved in their annihilation decays. We review this effective 
field theory and its application to heavy-quarkonium physics. Then we discuss the space- 
time structure of the QQ annihilation process and develop a general factorization formula 
for the annihilation decay rates of heavy quarkonia in terms of matrix elements of NRQCD. 

A. Energy Scales in Heavy Quarkonium 

In a meson containing a heavy quark and antiquark, there are several different momentum 
scales that play important roles in the dynamics. The most important scales are the mass M 
of the heavy quark, its typical 3- momentum Mv (in the meson rest frame), and its typical 
kinetic energy Mv"^. The heavy-quark mass M sets the overall scale of the rest energy of the 
bound state and also provides the short- distance scale for annihilation processes. The size 
of the bound state is the inverse of the momentum Mv, while Mv"^ is the scale of the energy 
splittings between radial excitations and between orbital-angular-momentum excitations. 
Spin splittings within a given radial and orbital-angular-momentum excitation are of order 
Mf^, but this scale plays no significant role in the dynamics. 

The typical velocity v of the heavy quark decreases as the mass M increases. If M is 
large enough, v is proportional to the running coupling constant ^^(M), and it therefore 
decreases asymptotically like 1/ log(M). Thus, if M is sufficiently large, the heavy quark and 
antiquark are nonrelativistic, with typical velocities v <^ 1. We assume in this paper that 
the mass M is heavy enough that the momentum scales M, Mv and Mv"^ are well-separated: 
{Mv'^Y ^ {MvY ^ M'^ ■ Quark potential model calculations indicate that the average value 
of f ^ is about 0.3 for charmonium and about 0.1 for bottomonium |[ll|, and these estimates 
are confirmed by lattice QCD simulations. Thus, the assumption {Mv'^Y *^ {MvY ^ M'^ 



is very good for bottomonium, and reasonably good even for charmonium. For lighter 
quarkonium states, such as the ss system, our analysis does not apply. 

Another momentum scale that plays a role in the physics of heavy quarkonium is Aqcd, 
the scale associated with nonperturbative effects involving gluons and light quarks. It de- 
termines, for example, the long-range behavior of the potential between the heavy quark 
and antiquark, which is approximately linear, with a coefficient of (450 MeV)^ |jTl|. We 
can use this coefficient as an estimate for the nonperturbative scale: Aqcd ~ 450 MeV. 
For both charmonium and bottomonium, the first radial excitation and the first orbital- 
angular-momentum excitation are both about 500 MeV above the ground state. Taking this 
value as an estimate for the scale Mv"^, we see that Aqcd and Mv"^ are comparable for both 
charmonium and bottomonium. 

Our analysis of heavy quarkonium annihilation is based on separating the effects at 
the momentum scale M from those at the lower momentum scales Mv, Mv"^ and Aqcd- 
The effects at the scale M are taken into account through the coupling constants of 4- 
fermion operators in the lagrangian for NRQCD. We assume that as{M) <^ 1, so that these 
coupling constants can be calculated using perturbation theory in as{M). The assumption 
that as{M) <^ 1 is well-satisfied for bottomonium, for which as{M) ^ 0.18, and reasonably 
well-satisfied for charmonium, for which as{M) ^ 0.24. 

The effects of the lower momentum scales Mv, Mv'^, and Aqcd are factored into matrix 
elements that can be calculated using nonperturbative methods, such as lattice-QCD simu- 
lations. These matrix elements are organized into a hierarchy in terms of their dependence 
on V. Our final expression for the annihilation rate therefore takes the form of a double 
expansion in as{M) and v. These expansion parameters are not independent for quarko- 
nium. The typical velocity v of the heavy quark is determined by a nonperturbative balance 
between its kinetic energy Mv'^/2 and the potential energy, which, for sufficiently large M, 
is dominated by a color-Coulomb term proportional to as{l/r)/r. Setting r ~ l/{Mv) in 
the potential and equating it with the kinetic energy, we obtain the identification 



V ~ as{Mv). (2.1) 

This equation can be solved self-consistently to obtain an approximate value for the typical 
velocity v. The identification ( p.l|) has a simple, but important, implication for calculations 
of annihilation rates. Since the running coupling constant in QCD decreases with the mo- 
mentum scale, V is greater than or of order as{M). Thus relativistic corrections of order 
(t;^)" can be expected to be more important than perturbative corrections of order a^"(M). 
In particular, there is little to be gained by calculating perturbative corrections at next-to- 
next-to-leading order in as{M), unless relativistic corrections through relative order v"^ are 
included as well. 

B. The NRQCD Lagrangian 

The most important energy scales for the structure and spectrum of a heavy quarkonium 
system are Mv and Mv"^, where M is the mass of the heavy quark Q and f <^ 1 is its average 
velocity in the meson rest frame. Momenta of order M play only a minor role in the complex 
binding dynamics of the system. We can take advantage of this fact in our analysis of heavy- 
quark mesons by modifying QCD in two steps. 

We start with full QCD, in which the heavy quarks are described by 4-component Dirac 
spinor fields. In the first step, we introduce an ultraviolet momentum cutoff A that is of 
order M. This cutoff explicitly excludes relativistic heavy quarks from the theory, as well as 
gluons and light quarks with momenta of order M. It is appropriate to our analysis of heavy 
quarkonium, since the important nonperturbative physics involves momenta of order Mv 
or less. Of course, the relativistic states we are discarding do have some effect on the low- 
energy physics of the theory. However, any interaction involving relativistic intermediate 
states is approximately local, since the intermediate states are necessarily highly virtual 
and so cannot propagate over long distances. Thus, generalizing standard renormalization 
procedures, we systematically compensate for the removal of relativistic states by adding 
new local interactions to the lagrangian. To leading order in 1/A or, equivalently, 1/M, these 
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new interactions are identical in form to interactions already present in the theory, and so 
the net effect is simply to shift bare masses and charges. Beyond leading order in 1/M, one 
must extend the lagrangian to include nonrenormalizable interactions that correct the low- 
energy dynamics order-by-order in 1/M. In this cutoff formulation of QCD, all effects that 
arise from relativistic states, and only these effects, are incorporated into renormalizations 
of the coupling constants of the extended lagrangian. Thus, in the cutoff theory, relativistic 
and nonrelativistic contributions are automatically separated. This separation is the basis 
for our analysis of the annihilation decays of heavy quarkonia. 

The utility of the cutoff theory is greatly enhanced if, as a second step, a Foldy- 
Wouthuysen-Tani transformation [|12| is used to block-diagonalize the Dirac theory so as 
to decouple the heavy quark and antiquark degrees of freedom. Such a decoupling of parti- 
cle and antiparticle is a familiar characteristic of nonrelativistic dynamics and is quite useful 
in our study of heavy quarkonium. The net effect is that the usual relativistic field theory of 
four-component Dirac spinor fields is replaced by a nonrelativistic Schrodinger field theory, 
with separate two-component Pauli spinor fields for the heavy quarks and for the heavy 
antiquarks. This field theory is NRQCD |^. The lagrangian for NRQCD is 

-^NRQCD = -flight + -^heavy + ^C,. (2.2) 

The gluons and the nf flavors of light quarks are described by the fully relativistic lagrangian 

Aight = -UiG.^G^'' + J^g^^g, (2.3) 

where G^i, is the gluon field-strength tensor expressed in the form of an SU(3) matrix, and q 
is the Dirac spinor field for a light quark. The gauge- covariant derivative is D^ = d'^ + igA'^, 
where A^ = (0, A) is the SU(3) matrix-valued gauge field and g is the QCD coupling 
constant. The sum in ( p.3| ) is over the nj flavors of light quarks. The heavy quarks and 
antiquarks are described by the term 

£...., = V* ('D, + ^) i. + x' U - ^) X, (2.4) 



where ip is the Pauh spinor field that annihilates a heavy quark, x is the Pauli spinor field 
that creates a heavy antiquark, and Dt and D are the time and space components of the 
gauge-covariant derivative D^. Color and spin indices on the fields ip and x have been 
suppressed. The lagrangian /^ught + -^heavy describes ordinary QCD coupled to a Schrodinger 
field theory for the heavy quarks and antiquarks. The relativistic effects of full QCD are 
reproduced through the correction term SC in the lagrangian (p.2|) . 

The correction terms in the effective lagrangian for NRQCD that are most important for 
heavy quarkonium are bilinear in the quark field or the antiquark field: 
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(V^t(D .gE-gE- B)^ + x^(D -gE-gE-B] 



X 



8M2 
+ ^ (V^^(^D xgE-gEx iB) ■ a^, + x^(^D x gE - gE x iB) ■ ax) 

+ ^ {^\g^ ■ ^)^ - x\g^ ■ ^)x) , (2.5) 

where E^ = G^* and B^ = ^e^^^G^^ are the electric and magnetic components of the gluon 
field strength tensor G^'^ . By charge conjugation symmetry, for every term in ( |2.5| ) involving 
if), there is a corresponding term involving the antiquark field x-, with the same coefficient q, 
up to a sign. The operators in ( [^.5D must be regularized, and they therefore depend on the 
ultraviolet cutoff or renormalization scale A of NRQCD. The coefficients Cj(A) also depend 
on A in such a way as to cancel the A-dependence of the operators. Renormalization theory 
tells us that NRQCD can be made to reproduce QCD results as accurately as desired by 
adding correction terms to the lagrangian like those in (|2.5|) and tuning the couplings to 



appropriate values 13 



Mixed 2-fermion operators involving x^ and ip (or ip'^ and x) correspond to the annihila- 
tion (or the creation) of a QQ pair. Such terms are excluded from the lagrangian as part of 
the definition of NRQCD. If such an operator annihilates a QQ pair, it would, by energy con- 
servation, have to create gluons (or light quarks) with energies of order M. The amplitude 
for annihilation of a QQ pair into such high energy gluons cannot be described accurately 
in a nonrelativistic theory such as NRQCD. Nevertheless, as is discussed in Section [II E| , the 
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effects of such annihilation processes on low energy amplitudes can be reproduced by adding 
4-ferniion operators such as '^/'^XX^V' to the effective lagrangian. 

Operators containing higher-order time derivatives, such as ip'^Dfip, are also omitted 
from the effective lagrangian as part of the definition of NRQCD. These operators can be 
eliminated by field redefinitions that vanish upon use of the equations of motion. Because 
of these field redefinitions, the off-shell Green's functions of NRQCD need not agree with 
those of full QCD, but the two theories are equivalent for on-shell physical quantities. 

The coefficients q in ( p.5|) must be tuned as functions of the coupling constant as, 
the heavy-quark mass parameter in full QCD, and the ultraviolet cutoff A of NRQCD, so 
that physical observables are the same as in full QCD. The coefficients are conveniently 
determined by matching low-energy scattering amplitudes of heavy quarks and antiquarks 
in NRQCD, calculated in perturbation theory in a^ and to a given precision in v, with the 
corresponding perturbative scattering amplitudes in full QCD. It is necessary to use on-shell 
scattering amplitudes for this purpose, because the equations of motion have been used to 
simplify the effective lagrangian for NRQCD by eliminating terms with more than one power 
of Dt. The scattering amplitudes can be calculated using perturbation theory in a^, since 
the radiative corrections to the coefficients in the NRQCD lagrangian are dominated by 
relativistic momenta. These coefficients therefore have perturbative expansions in powers of 
as{M) P,p!4|. The coefficients in ( |2.5|) are defined so that Cj = 1 + 0{as). 

The explicit factors of M in (|2.5|) were introduced in order that the coefficients c, be 
dimensionless. These coefficients therefore depend on the definition of the heavy-quark 
mass parameter M. Our definition of M is specified by the lagrangian (|2.4| ): 1/(2M) is the 
coefficient of the operator ip'^'D'^ip. If a different prescription is adopted for M, then all the 
Cj's must be changed accordingly. The simplest way to determine the mass parameter M is to 
match the location of the pole in the perturbative propagator for a heavy quark in NRQCD 
with that in full QCD. In both NRQCD and full QCD, the kinetic energy for a heavy quark of 
momentum p in perturbation theory has the form E = p^/(2Mpoio) —p'^/{8M^^iJ + . . ., where 
Mpoie is the perturbative pole mass. In Appendix p3 1| , the self-energy of the heavy quark 
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is calculated in NRQCD to order a^ and to leading order in v. If we use a regularization 
scheme in which power divergences are subtracted, then the energy-momentum relation gives 
M = Mpoie(H-0(af)). The corresponding calculation using a lattice regularization has been 
carried out by Morningstar [|T5[. The perturbative pole mass can be related to any other 
definition of the heavy-quark mass by a calculation in full QCD. 

C. Velocity-scaling Rules 

In principle, infinitely many terms are required in the NRQCD lagrangian in order to 
reproduce full QCD, but in practice only a finite number of these is needed for precision 
to any given order in the typical heavy-quark velocity v. We can assess the relative im- 



portance of various terms by using velocity-scaling rules that were derived in Ref. |TJ] and 
are summarized in Table |. This table lists the fields and operators from which terms in 
the NRQCD action are built, together with the approximate magnitude of each for matrix 
elements between heavy quarkonium states that are localized in space. The scaling rules 
were derived in Ref. []T^ by analyzing the equations of motion for the quantum field opera- 
tors of NRQCD. The typical heavy-quark velocity v is determined dynamically by a balance 
between the kinetic and potential terms in the equation of motion for the heavy-quark field, 
and V can be used as an expansion parameter in order to analyze the importance of other 
terms. The scaling rules are certainly correct within perturbation theory in a^, but, since 
they are based on the self-consistency of the field equations, they should also be valid in the 
presence of nonperturbative effects. 

There is an important caveat to the velocity-scaling rules that involves ultraviolet- 
divergent loop corrections. Loop corrections to an operator give rise to power ultravio- 
let divergences, as well as to logarithmic divergences. The logarithmic divergences modify 
the scaling rules by factors of log(A/Mf). The power divergences can contribute factors of 
1/f ", and the scaling rules apply only after such 1/f " divergences have been subtracted. The 
subtracted expression is the relevant one for the following reason. The power-divergent con- 
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tributions to a given operator O that yield factors of 1/f " have the form of renormahzations 
of lower-dimension operators. When the coefficients of NRQCD are tuned so as to repro- 
duce full QCD, the coefficients of the lower- dimension operators are adjusted so that their 
contributions to physical quantities cancel the contributions of the 1/f" power-divergent 
loop corrections to the operator O. Consequently, the inclusion of a given operator in the 
NRQCD lagrangian yields a net correction to any physical quantity that is in accordance 
with the velocity-scaling rules, up to logarithmic corrections. 

The estimates for the magnitudes of g(j) and gA in Table | hold in Coulomb gauge. 
Coulomb gauge is a natural gauge for analyzing heavy quarkonium, because it avoids spu- 
rious retardation effects that are present in covariant gauges, but cancel out in physical 
quantities [|1^. Coulomb gauge is also a physical gauge, that is, a gauge with no negative 
norm states. Thus, it allows a sensible Fock-state expansion for the meson. The dominant 
Fock state is of course \QQ), but the meson also contains the Fock state \QQg), which 
includes a dynamical gluon, and higher Fock states as well. 

The estimates in Table | were derived assuming that one can do perturbation theory in 
the typical heavy-quark velocity. This perturbation theory relies on the fact that soft gluons 
have a weak coupling to heavy quarks, not because the coupling constant as is small, but 
because the interaction is proportional to the heavy-quark velocity v. In the derivation of 
the magnitude of gA in Ref. |]14|, dynamical gluons were assumed to have typical momenta 



of order Mv, which is the inverse size of the quarkonium. The perturbative estimate for the 
magnitude of the operator gA is as{k)vk for a dynamical gluon of momentum k. For k of 
order Mv, we can set as ^ v and recover the estimate Mv^ given in Table |. For k of order 
Mv^, we can set a^ ~ 1, and we again obtain the estimate Mv^. This estimate relies on 
perturbation theory, which may be suspect because of the strong coupling between gluons 
with momenta on the order of Mv'^. However such gluons necessarily have wavelengths of 
order l/[Mv'^) or larger, which is much larger than the typical size l/{Mv) of the quarko- 
nium. For such long-wavelength gluons, the multipole expansion, whose validity transcends 
that of perturbation theory in the coupling constant, can be used to justify the estimate for 
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^A in Table 10. 

The velocity-scaling rules in Table | show that the terms in SCuimear in ( P-5| ) all give con- 
tributions that are suppressed by 0{v^) relative to those from the leading lagrangian /^heavy- 
Recalling that mixed 2-fermion operators, such as ?/'^(D^)^x, and operators involving higher 
time derivatives, such as ip^D^ip, are omitted as part of the definition of NRQCD, we see that 
SCuiincai contains all the 2-fermion NRQCD operators of relative order v"^. The lagrangian 
/flight + 'Cheavy + (^ /^bilinear cau therefore bc used to calculate NRQCD matrix elements be- 
tween heavy-quarkonium states with an error of order v^. If an error of order f ^ is sufficiently 
accurate, then the matrix elements can be calculated by using the lagrangian ^ught + /^heavy- 
It is instructive to contrast the relative magnitudes of the NRQCD operators in the case 
of a heavy quarkonium with the relative magnitudes of the same operators in the case of 
a heavy-light meson. (In the meson rest frame, the lagrangian for NRQCD is identical to 
that for Heavy Quark Effective Theory, which is the standard formalism for treating heavy- 
light mesons [1^.) In a heavy-light meson, the typical 3-momentum of the heavy quark is 
of order Aqcd, and is independent of the heavy-quark mass. The binding energy is also 
of order Aqcd, and is much larger than the heavy-quark kinetic energy, which is of order 
Aqcd/M. Thus, in a heavy-light meson, the 3-momentum and the energy of the heavy-quark 
are both of order Aqcd, in contrast with the situation in a heavy quarkonium, in which the 
3-momentum is of order Mv and the energy is of order Mv'^. Consequently, in a heavy-light 
meson, the effects of operators of dimension d are of order {Aqcd/^Y"'^ relative to the 
effects of the dimension-4 operator ipHDtip- The leading term ipHDtip describes a static 
heavy quark acting as a source of gluon fields. All effects of relative order Aqcd/M can be 
taken into account by adding the dimension-5 operators ip'^'D'^il! and ip'^gB ■ crip. 

D. Quarkonium in NRQCD 

Several qualitative features of heavy quarkonium can be inferred directly from the 
NRQCD lagrangian by exploiting the heavy-quark velocity v as an expansion parameter. 
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Expansions in powers of v are possible in ordinary QCD, but they are complicated by the 
need to make a nonrelativistic expansion of each individual Lorentz-invariant operator in 
order to separate the various powers of v. Relativistic effects have been unraveled to a large 
extent in NRQCD, with the leading f-dependence of each operator being specified by the 
velocity-scaling rules in Table |. 

The most distinctive phenomenological feature of heavy quarkonium is that, for many 
purposes, it is accurately described by the quark potential model, in which the heavy quark 
and antiquark are bound by an instantaneous potential. This model is a tuned phenomenol- 
ogy, rather than a theory, but it is far simpler than a full field-theoretic description based 
on NRQCD or QCD. Its validity rests upon two essential ingredients of heavy- quarkonium 
physics. The first is that the dominant effect of the exchange of gluons between the heavy 
quark and antiquark is to produce an instantaneous interaction. The reason for this is that 
the most important gluons have momenta of order Mv and energies of order Mv"^. Such 
gluons are off their energy shells by amounts of order Mv, which are much greater than the 
typical kinetic energy Mv'^ of the heavy quark. Consequently, the interaction times of the 
gluons are shorter by a factor of 1/v than the time scale associated with the motion of the 
heavy quarks, and the gluons' interactions are, therefore, instantaneous as far as the heavy 
quarks are concerned. 

The second essential ingredient underlying the quark potential model is that the prob- 
ability of finding dynamical gluons (those that are not part of the potential) in the meson 
is small. This is important because dynamical gluons with very low energy produce effects 
that are not instantaneous and are not readily incorporated into the quark potential model. 
In particular, gluons with energies of order Mv'^ have interaction times comparable to that 
of the heavy quarks, and their exchange therefore leads to significant retardation effects. 
The probability for the Fock state \QQg) of the meson can be estimated by considering the 
energy shift of a quarkonium state \H) that is due to the presence of a Fock-state component 
\QQg)- In Coulomb gauge, the only terms that connect the dominant Fock state \QQ) to 
the Fock state \QQg) are terms that involve the vector potential A. At leading order in v, 
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the contributions to the energy shift come from the term igA ■ ip^'Vip/M in £heavy: 

AE = -ITT {H\ f d^xigA-ij^Vij\H). (2.6) 

Using the velocity-scahng rules in Table |I| and taking into account the relevant integration 
volume l/{Mv)^, we obtain the estimate AE ~ Mf^. This energy shift can be written 
in a different way — as the product of the probability -Pqq„ for the QQg state multiplied 
by the energy Eqq of that state. For gluons with momenta k of order Mv, the energy 
of the QQg state is dominated by the energy of the gluon, and we find that PQQg ~ v^. 
For dynamical gluons with very low energies of order Mv"^ or less, the energy of the QQg 
state is of order Mt>^ and we obtain the estimate Pnog ~ v^. For heavy quarkonium, QQg 
states are therefore suppressed relative to the dominant QQ state by a factor of order f ^ 
in the probability. Hence, for most quantities, effects due to Fock states like \QQg) that 
contain dynamical gluons are suppressed by powers of v. This might be expected from the 
phenomenological successes of the quark potential model. However, there are quantities, 
such as the decay rates of P-wave states into light hadrons |T^, for which the effects of the 
Fock state \QQg) are of leading order in v and the quark potential model fails completely. 

The above estimates for the probabilities of \QQg) Fock states apply if the spin state of 
the QQ pair is the same as in the dominant \QQ) Fock state. If the spin state is different, 
we must replace gA ■ V in (2.6) with (7B ■ a to obtain a nonzero matrix element. Using 
the velocity-scaling rules of Table I, we again obtain an estimate AE ~ Mv^ for the energy 
shift, implying that the probability for a \QQg) state containing a gluon with momentum 
on the order of Mv is -Pgog ~ v^- However, in the derivation of the velocity-scaling rules in 
Ref. [14], it was assumed that dynamical gluons have momenta of order Mv. If the gluon has 
a much smaller momentum fc, then the estimate M^f ^ for the operator (7B in Table I should 
be replaced with k'^v'^. Using this to estimate the energy shift from a \QQg) Fock state 
containing a gluon with momentum of order Mf^, we obtain AE ~ Mv^ and Pn-Qq ~ v^- 
Thus, gluons with very low momenta exhibit the suppression that is characteristic of the 
multipole expansion. We conclude that a \QQg) Fock state that can be reached from the 
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dominant \QQ) Fock state by a spin- flip transition is dominated by dynamical gluons with 
momenta of order Mv and that the probability of such a Fock state is Pqq„ ~ v^. 

Another important feature of quarkonium structure is its approximate independence 
of the heavy-quark spin. This feature follows immediately from the structure of the 
NRQCD lagrangian, which exhibits an approximate heavy-quark spin symmetry. The lead- 
ing term £heavy is completely independent of the heavy-quark spin. With just this term, 
states that differ only in the spins of the heavy quark and antiquark have identical proper- 
ties; heavy-quark spin is conserved and can be used to label the energy eigenstates. Spin- 
dependence enters first through the bilinear terms in (|2.5| ) that contain Pauli matrices, and 
they give corrections that are of relative order v^B Thus, spin splittings for quarkonia should 
be smaller than splittings between radial and orbital-angular-momentum excitations, with 
the ratios of these splittings scaling roughly as f^. This familiar feature of the spectra of 
charmonium and bottomonium reinforces our confidence in the power-counting rules and in 
the utility of a nonrelativistic framework for studying quarkonium. 

The total angular momentum J, the parity P, and the charge conjugation C are exactly 
conserved quantum numbers in NRQCD, as well as in full QCD. Thus, the energy eigenstates 
\H) of heavy quarkonium can be labelled by the quantum numbers J^*^. By the arguments 
given above, the dominant component in the Fock state expansion of \H) is a pure quark- 
antiquark state \QQ). The Fock state \QQg), in which a dynamical gluon is present, has a 
probability of order v"^, and higher Fock states have probabilities of order f ^ or higher. Since 
our primary interest is in processes in which the Q and Q in the quarkonium annihilate. 



^In perturbation theory, ladder-like Coulomb-gluon exchanges between the quark and antiquark 
give a factor of order Ug/v for each ladder rung. The spin- flip contribution is down by f ^ relative 
to this Coulomb-ladder contribution. For example, in a two-loop calculation, the Coulomb ladder 
gives a factor of order {os/v)'^, while the ladder with one Coulomb exchange and one spin-flip 
exchange gives a factor of order v'^{as/v)'^ = a^. 
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we concentrate on the state of the QQ pair in the various Fock-state components. For a 
general Fock state, the QQ pair can be in either a color-singlet state or a color-octet state. 
Its angular- momentum state can be denoted by the spectroscopic notation '^'^^^Lj, where 
S* = 0, 1 is the total spin of the quark and antiquark, L = 0, 1, 2, . . . (or L = S,P,D,...) 
is the orbital angular momentum, and J is the total angular momentum. A QQ pair in a 
"^^^^Lj state has parity P = (—1)^^"'^; if it is in a color-singlet state, it has charge-conjugation 
number C = (-1)^+'^. 

In the Fock state \QQ), the QQ pair must be in a color-singlet state and in an angular- 
momentum state "^^^^Lj that is consistent with the quantum numbers J^'" of the meson. 
Conservation of J^'" implies that mixing is allowed only between the angular-momentum 
states ^(J — l)j and ^(J + l)j. For example, a ^5*1 QQ state can mix with a ^Di state. 
However, such mixing is suppressed because operators that change the orbital angular mo- 
mentum must contain at least one power of V. In general, up to corrections of order f^, 
we can regard the QQ component of the meson as being in a definite angular-momentum 
state ^^^^Lj. Of course, if the contribution of the dominant angular-momentum state is 
suppressed in a given process, then the contribution of the subdominant states takes on 
increased importance. We will present examples of this phenomenon in the discussions of 
the decay and production of P-wave states. 

We turn next to the Fock state \QQg) of the meson, which includes a dynamical gluon 
and has a component whose probability is of order v'^. In spite of the fact that the dynamics 
of the soft gluon is nontrivial, NRQCD tells us much about the quantum numbers of the 
QQ pair in the QQg component whose probability is of order v"^. The pair must of course 
be in a color-octet state. Heavy-quark spin symmetry implies that the total spin quantum 
number S for the QQ pair is the same as in the dominant Fock state \QQ). But NRQCD 
also tells us that the orbital state of the QQ pair is closely related to that in the Fock 
state \QQ). The reason for this is that the coupling of the soft gluon can be analyzed 
using a multipole expansion, and the usual selection rules for multipole expansions apply. 



The leading interaction that couples the dominant Fock state \QQ) to the state \QQg) is 
the electric-dipole part of the operator ■ip'^gA ■ V-?/' in £heavy, and this changes the orbital- 
angular-momentum quantum number L of the QQ pair by ±1. Higher multipoles bring in 
additional powers of f , as does second-order perturbation theory. Thus, if the QQ pair in 
the dominant Fock state \QQ) has angular-momentum quantum numbers '^^^^Lj, then the 
Fock state \QQg) has a probability of order v"^ only if the QQ pair has total spin S and 
orbital angular momentum L + 1 or L — 1. For example, if the dominant Fock state consists 
of a QQ pair in a ^5*1 state, then the Fock state \QQg) has a probability of order v"^ only if 
the QQ pair is in a color-octet state with angular-momentum quantum numbers ^Pq, "^Pi, 
or ^P2- If the dominant Fock state consists of a QQ pair in a ^Pi state, then the Fock state 
\QQg) has a probability of order v"^ only if the QQ pair is in a color-octet ^5*0 or ^D2 state. 
The above discussion applies to Fock states \QQg) in which the QQ pair has the same 
total spin quantum number S as in the dominant \QQ) state. The probabilities for Fock 
states \QQg) that can be reached from the dominant Fock state by a spin-flip transition also 
scale in a definite way with v. The probability for such a Fock state to contain a dynamical 
gluon with momentum of order Mv is of order f ^, just as in the case of a non-spin-fiip 
transition. However, in the case of a spin-flip transition, this momentum region dominates 
because, as we have seen, gluons with softer momenta, on the order of Mf ^, are suppressed 
by the multipole expansion. Thus, if the QQ pair in the dominant Fock state has angular- 
momentum quantum numbers "^^^^Lj, then the Fock state \QQg), with the QQ pair in a 
color-octet state with the same value of L but different total spin quantum number, has a 
probability of order v^. For example, if the dominant Fock state consists of a QQ pair in 
a ^5*1 state, then the Fock state \QQg) with the QQ pair in a color-octet ^Sq state has a 
probability of order v^. If the dominant Fock state consists of a QQ pair in a ^ Pi state, then 
the Fock state \QQg) with the QQ pair in a color-octet ^Pj state has probability of order 



t>^. 



19 



E. Space-Time Structure of Annihilation 

As we will explain in this subsection, the annihilation of a heavy QQ pair into gluons 
(or light quarks) occurs at distances that are typically of order 1/M, that is, at momentum 
scales of order M. Because of the large momentum scales involved, the details of the annihi- 
lation process cannot be described accurately within a nonrelativistic effective theory such 
as NRQCD. Nevertheless, as we will argue in the next subsection, the effects of annihila- 
tion can be incorporated into NRQCD through 4-fermion operators in the term 5C in the 
NRQCD lagrangian. To show that the required operators are local, it is sufficient to show 
that the interactions they account for occur over short distances of order 1/M. Strictly local 
operators are then obtained by expanding the short-distance interaction in a Taylor series 
in the 3- momentum p of the heavy quark multiplied by the characteristic size 1/M. 

Now we wish to argue that the annihilation process is indeed local, i.e. that the anni- 
hilation does occur within a distance of order 1/M. We note that any annihilation must 
result in at least two hard gluons (or light quarks), each with momentum of order M. This 
has two consequences. First, the heavy quark and antiquark must come within a distance 
of order 1/M in order to annihilate. That is because the emission of a hard gluon from, 
say, the heavy quark puts it into a highly virtual state, which can propagate only a short 
distance before the quark must annihilate with the antiquark. Thus, the total annihilation 
amplitude can be expressed as the sum of point-like annihilation amplitudes, where the sum 
extends over the possible annihilation points inside the meson. The annihilation rate is 
the square of the total annihilation amplitude, summed over all possible final states. The 
second consequence of the hard gluons is that there is no overlap between one annihilation 
amplitude and the complex conjugate of another if the two annihilation points are sepa- 
rated by a distance greater than about 1/M. This might seem surprising, since the gluons 
are, in effect, on their mass shells (that is, they fragment into jets with invariant masses 
much less than M). There is no highly virtual state to constrain the distance between the 
annihilation points for two amplitudes that produce the same final-state jets. Nevertheless, 
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the annihilation points must be in close proximity to each other in order for there to be an 
overlap between the final states. In order to see why this is so, we note that, in classical 
mechanics, we could trace the two final-state jets back to the annihilation vertex, and there 
would be no ambiguity whatsoever as to its space-time position. In quantum mechanics, 
the uncertainty principle tells us that we can know the position of the annihilation vertex 
only to a precision of order 1/M, since the jet momenta are of order M . Hence, in quantum 
mechanics, QQ annihilation is not a point-like process, but it is a localized process, with a 
size of order 1/M. 

In a field-theoretic calculation of the annihilation rate at leading order in a^, the local- 
ization of the annihilation process would manifest itself as follows. The annihilation rate 
involves the imaginary part r{P,p,p') of the scattering amplitude for a QQ pair with to- 
tal momentum P and initial and final relative momenta p and p'. Consider the Fourier 
transform of r{P,p,p') with respect to all three momentum variables: 

j (fPd^pdy gip.(x-x')gip-(xi-.2)g^p'.K-4) Y{P,p,p'). (2.7) 

Here, Xi and X2 correspond to quark and antiquark interaction points in one annihilation 
amplitude, x'^ and x'2 correspond to quark and antiquark interaction points in the complex 
conjugate of a second annihilation amplitude, and X = {xi + X2)/2 and X' = {x'l + x'^)!'^ 
are average annihilation points for the first and second amplitudes. The fact that r{P,p,p') 
is insensitive to changes in p and p' that are much less than M implies that, in the Fourier 
transform, Xi {x[) is localized to within a distance of order 1/M of X2 (3^2)- Similarly, the 
fact that r is insensitive to changes in P that are much less than M implies that the first and 
second amplitudes have significant overlap only if X and X' are separated by a distance of 
order 1/M or less. Note that, if one puts a restriction in the annihilation rate on one of the 
components of the jet momentum, then F becomes sensitive to that component of P, and 
the annihilation vertices are no longer localized along that direction. This is a consequence 
of the uncertainty principle, which says that knowledge of a component of the jet momentum 
along a given direction reduces our potential knowledge of the position of the annihilation 
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vertex along that direction. 

The radiation of soft or colhnear gluons might seem to violate this simple localization 
picture that appears at leading order in the coupling constant. Gluon radiation from the 
initial QQ pair is not a problem, since infrared divergences can be factored into the long- 
distance matrix elements of the 4-fermion operators that mediate the annihilation process 
in NRQCD, and collinear divergences are controlled by the heavy-quark mass. We must, 
however, worry about infrared or collinear divergences from the radiation of gluons from 
the final-state hard gluons. In the presence of such soft or collinear radiation, the hard 
gluon can propagate almost on its mass shell from the annihilation point to the emission 
vertex. The energetic final-state gluon jet points back to the emission vertex, rather than 
to the annihilation point, which may be far away. In perturbation theory, infrared and 
collinear divergences occur in individual Feynman diagrams and produce a sensitivity to 
the heavy-quark momenta in F. However, the Kinoshita-Lee-Nauenberg (KLN) theorem 
p!9| guarantees that, when one sums over the contributions of all nearly degenerate final 
states, as is done in forming the inclusive annihilation rate, the infrared divergences cancel 
between diagrams involving real and virtual gluon emission. We can think of this KLN 
cancellation as a consequence of a generalized form of the uncertainty principle: we can 
localize the annihilation point, provided that we do not require too much knowledge about 
the final state — that is, provided that we do not distinguish between the various states that 
contribute to the inclusive cross section. 

The locality of the annihilation process is spoiled if the final-state gluons form a narrow 
resonance, such as a glueball. This is because the jets produced by the decay of the resonance 
point back to the place where the resonance decayed. If the resonance is narrow, this 
may be far from the point where the heavy quark and antiquark annihilated. That is why 
perturbation theory cannot be applied directly to the cross section for e+e~ annihilation into 
hadrons in the region of the charmonium or bottomonium resonances. In a field-theoretic 
calculation, the resonance partially spoils the KLN cancellation of infrared and collinear 
divergences. While contributions from gluons that have exactly zero momentum or are 
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exactly collinear still cancel, the real and virtual contributions no longer cancel for soft gluons 
whose energy is comparable to the resonance width or collinear gluons whose transverse 
momentum is comparable to the resonance width. In the case of e^e" annihilation, one can 
deal with this problem by forming a suitable average of the cross section over the resonance 
region |2^. In perturbation theory, the effect of this smearing is to allow virtual soft or 



collinear emission at one value of the e+e^ center-of-mass energy y/s to cancel real soft or 
collinear emission at a slightly higher value of ^/s, but the same value of the energy of the 
resonating QQ pair. This solution of smearing in the energy is not available to us in the case 
of quarkonium annihilation. Fortunately, there are no known narrow glueball resonances in 
the charmonium or bottomonium region, so we do not expect the resonance issue to be a 
problem in practice. 

F. Annihilation into Light Hadrons 

Since the annihilation of a QQ pair necessarily produces gluons or light quarks with ener- 
gies of order M , the annihilation amplitude cannot be described accurately within NRQCD. 
Nevertheless, the annihilation rate, which is the square of the amplitude summed over fi- 
nal states, can be accounted for in NRQCD. Since the annihilation rate of the QQ pair 
is localized within a distance of order 1/M, the annihilation contribution to a low-energy 
QQ —>■ QQ scattering amplitude can be reproduced in NRQCD by local 4-fermion operators 
in 6C involving ip, x\ Xi ^"^^ V"^- The optical theorem relates QQ annihilation rates to 
the imaginary parts of QQ -^ QQ scattering amplitudes. This relation implies that the 
coefficients of the 4-fermion operators in 5C must have imaginary parts. These imaginary 
parts are the manifestation of annihilation in NRQCD. 

The 4-fermion interactions that represent the effects of QQ annihilation in NRQCD have 
the general form 

^A-fermion = ^ »L -4 ^"(A), (2.8) 

n -'" 
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where the On are local 4-fermion operators, such as ip^xx^i'- The naive scaling dimensions 
dn of the operators can be obtained by counting the powers of M using Table |. The factors 
of M'^"~^ in ( |2.8D have been introduced so as to make the coefficients /„ dimensionless. The 
operators On must be regularized, and they therefore depend on the ultraviolet cutoff or 
renormalization scale A of the effective theory. The natural scale for this cutoff is M, since 
1/M is the distance scale of the annihilation process. However, all results are independent 
of A, since the coefficients /n(A) depend on A in such a way as to cancel the A-dependence 
of the operators. The coefficients can be computed in full QCD as perturbation series in 
as{M), in which individual terms may depend on log(M/A). 

If the analysis of annihilation rates were carried out completely within full QCD, then 
the scale A would arise as an arbitrary factorization scale that must be introduced in order 
to separate the momentum scale M from smaller momentum scales of order Mv or less. 
The factorization scale A should not be confused with the renormalization scale /x of the full 
theory. The coefficients /„(A) are independent of fi if they are computed to all orders in 
as(/i), although some //-dependence is introduced as usual by the truncation of the pertur- 
bation series. Unless we explicitly specify otherwise, we always make the choice // = M in 
this paper. 

The dimension-6 4-fermion terms in 6C are 

+ ^^OsCS,) + ^OsCS,), (2.9) 

where the dimension-6 operators are 

OiCSo) = i^hx^i^, (2.10a) 

Oi(35i) = ^Vx-X^o-^, (2.10b) 

08('5o) = V^^T»XX^T>, (2.10c) 

OsCSi) = ^^aT^x ■ xVT>. (2.10d) 

The subscript 1 or 8 on the operators and on their coefficients indicates the color structure of 
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the operator. The arguments "^^^^Lj indicate the angular-momentum state of the QQ pair 
which is annihilated or created by the operator. Normal-ordering of the 4-fermion operators 
On will always be understood, so that matrix elements of On receive contributions only from 
annihilation of the Q and Q. The dimension-8 terms in the lagrangian for NRQCD include 



{5C. 



4— fermionJ(^=8 



/i( -Pi) n (^p\ -u ^^^ ^°) n i^p\ -u ■^i^ -^i) n (^p\ 



fli'P2 



OA'P2) 



9ii^So) ^ ^le N , 9ii^Si) ^ 3 
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Vii'S,) 



(2.11) 



The dimension-8 operators included explicitly in ( p.ll| ) are 



Oi('Po) = l 



2-7A A V 2^)^^, 



Oi(3Pi) = - V^t(_iD X ct)x ■ X^(-fD X ct)^, 



O1CP2) 



^t(_iBv))^^t(_iB(v))^, 



\lj^Xx\-0?^ + h-c. 






^Vx ■ x^o-(-f D)2^ + h.c. 



(2.12a) 
(2.12b) 
(2.12c) 
(2.12d) 
(2.12e) 
(2.12f) 
(2.12g) 



where D is the difference between the covariant derivative acting on the spinor to the right 
and on the spinor to the left: x^D-?/; = x^(Di^) ~ (Dx)^V^- We have used the notation T^^^^ 
for the symmetric traceless component of a tensor: T^^^^ = (T^^ + T-'*)/2 — T^''5^^ /3. For 
each of the operators shown explicitly in ( |2.11|) , there is a corresponding color-octet operator 
Os or Vs, which contains color matrices T" inserted between ^/^^ and x ^i-nd between x^ and 
ip. This exhausts the list of the dimension-8 operators that contribute at tree level to QQ 
scattering in the center of momentum frame. There are other dimension-8 operators, such 
as V('?/'^x) ■ (x^D-?/;) or T){ip'^T"'x) ' D(x^7'°'?/'), in which a derivative acts on the product 
of ip"^ and X o^^ on the product of x^ and ip. Matrix elements of operators such as these 
are proportional to the total momentum of the QQ pair, and therefore do not receive any 
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contributions from the dominant Fock state \QQ) in the meson rest frame. They do, however, 
receive contributions from higher Fock states, such as \QQg), in which the total momentum 
of the QQ pair is nonzero. 

According to the velocity-scahng rules in Table |I|, the dimension-6 terms in ( |2.9| ) scale 
as V relative to the leading term £heavy in the NRQCD Lagrangian. Thus, if we consider 
only the dependence on v, the terms in ( p.9[ ) appear to be more important than the terms in 
(5£biiinear, which scalc as v^ relative to the terms in £heavy However, the contributions from 
4-fermion operators contain extra suppression factors, owing to the operator coefficients, 
whose imaginary parts are of order a1{M) or smaller. Thus, the contributions to annihilation 
widths from ( p.9|) are of order af (M)f or smaller relative to the scale Mf ^ of the splittings 
between radial excitations and between orbital-angular-momentum excitations. Similarly, 
the contributions to annihilation widths from the dimension-8 operators in (|2.11| ) are at most 



of order a'^{M)v^ relative to the scale Mv'^ of splittings between energy levels. Thus, the 
annihilation decay rates for heavy-quarkonium states are tiny perturbations on the energy 
levels. This is certainly true empirically. In the charmonium system, the ground state rjc has 
the largest annihilation width, but it is less than 3% of the splitting between the rjc and the 
first radial or orbital-angular-momentum excitations. For bottomonium, the annihilation 
widths are always less than 1% of the corresponding splittings. 

In order to obtain an expression for the annihilation rate, we recall that the decay 
rate is —2 times the imaginary part of the energy of the state. The contribution to the 
imaginary part of the energy that corresponds to annihilation into light hadrons comes from 
the expectation value of —^C^^fcrmion, whose coefficients have imaginary parts. Thus, we see 
that the annihilation rate of a heavy-quarkonium state H into light hadrons is 

r{H ^ LH) = 2 Im (i^l^A-fermionli^), (2.13) 

where LH represents all possible light-hadronic final states. The expectation value is taken 
in the rest frame of the quarkonium, where its total momentum P vanishes. The state \H) = 
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\H(P = 0)) is an eigenstate of the NRQCD hamiltonian.S It has the standard nonrelativistic 
normahzation: {H{P')\H{P)) = {27ry'^5^{P - P'). Inserting the expansion (^) into ( pl3D , 
we obtain 

r{H - LH) = E ^^5^ {H\OniA)\H). (2.14) 



The equation (|2.14| ) is our central result for the annihilation decays into light hadrons. It 
expresses the decay rate as a sum of terms, each of which factors into a short-distance coef- 
ficient Im/„ and a long-distance matrix element {H\On\H). The coefficients Im/„ in (p. 14 ) 



are proportional to the rates for on-shell heavy quarks and antiquaries to annihilate from 
appropriate initial configurations into hard gluons and light quarks, and can be computed 
as perturbation series in as{M). The matrix elements {H\On\H) give the probability for 
finding the heavy quark and antiquark in a configuration within the meson that is suitable 
for annihilation, and can be evaluated nonperturbatively using, for example, lattice simula- 
tions. The dependence on the arbitrary factorization scale A in ( 2.14 ) cancels between the 
coefficients and the operators. 



^Radial and orbital-angular-momentum excitations of a quarkonium may decay through the her- 
mitian part of the NRQCD lagrangian to lower-lying quarkonium states plus light hadrons. An 
example is the decay of ip{2S) into ipTrn. In this example, the spectrum of states in NRQCD con- 
tains a continuum of ^tttt scattering states, each of which includes a small admixture of the bare 
■0(25) state, and a discrete state, which is mostly the bare ip{2S) state, but which also contains a 
small admixture of bare ■07r7r scattering states. The ip{2S) Breit-Wigner resonance in, for example, 
the amplitude for e"'"e^ -^ A*^/^^ results from the contributions of the complete spectrum of states. 
However, the resonance in the amplitude can be reproduced by a single state, with complex energy, 
that is an eigenstate of the nonlocal effective Hamiltonian that one would obtain by integrating 
out the light-hadron states in NRQCD. One should identify the state \H) in ( |2.14 ) with such an 



eigenstate in applying (p. 14) to an excited quarkonium state that decays through the hermitian 



part of the NRQCD lagrangian into a lower-lying quarkonium state. 
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In some calculations of the matrix elements in ( p.l4|) , such as lattice simulations, it may 



be useful to approximate the states \H) by eigenstates of the hermitian part of the NRQCD 
Hamiltonian. We note that corrections to this approximation first appear at third-order in 
perturbation theory in Im(5£4_fcrmion, since second-order perturbation theory does not give 
an imaginary contribution to the energy. The corrections are therefore of order {T / Mv'^YV . 
This is of relative order a^(M)t'^ or smaller, since the leading terms in F scale like Mv^ 
and are multiplied by short-distance coefficients of order a1{M) or smaller. This level of 
accuracy is sufficient for most practical purposes. 

Applying the velocity-scaling rules of Table | to the matrix elements (if |(9„(A)|if ), one 
finds that the expression ( |2.14|) for the annihilation decay rate can be organized into an 
expansion in powers of v. Only a finite number of operators contribute to any given order 
in V. The coefficients fn{^) can be calculated as perturbation series in as{M), so (p.l4|) is 



really a double expansion in as{M) and v. The simultaneous expansion in as{M) and v is 
useful to the extent that these two parameters are both small. Of course, as{M) and v are 
not independent for heavy quarkonium. According to ( ^.11 ), v can be identified with as{Mv), 
which is larger than as{M). This implies that it would be futile to consider corrections to 
the coefficients Im/„ of relative order a'^{M) unless matrix elements {H\On\H) of relative 
order f " have already been included. 

The relation between v and as{M) implied by (|2.1|) follows from the dynamics of heavy 
quarkonium. The factorization formula (|2.14|) is actually an operator equation, and it can 
equally well be applied to other problems in which the relation between v and as{M) is 
different. An example in which v and as{M) are independent is the annihilation of a pair 
of heavy-light mesons, such as D and D mesons, at small relative velocity f <^ 1. As long 
as V is much larger than Kqcd/^Ii which is the typical relative velocity of a heavy quark 
in the heavy-light meson, it can be identified with the velocity of the heavy quark and the 
scahng rules of Table | apply. 
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G. Electromagnetic Annihilation 

In addition to annihilating into light hadrons, heavy-quarkonium states can also annihi- 
late into purely electromagnetic final states containing only photons and lepton pairs. The 
energies of the final-state photons and leptons are of order M. In NRQCD, the effects of 
electromagnetic annihilation can be accounted for in the same way as the effects of annihi- 
lation into light hadrons: by adding 4-fermion terms 5C^^^^^^:^„^ to the effective lagrangian. 
The primary difference is that in the case of electromagnetic annihilation, the final state, as 
far as the strong interactions are concerned, is the QCD vacuum state |0). The 4-fermion 
operators that reproduce the effects of electromagnetic annihilation therefore differ from 
those in (|2.9|) and ( p.ll|) by the insertion of an operator |0)(0| that projects onto the QCD 
vacuum state. The dimension-6 terms that must be added to the lagrangian are 

K"fe™),^, = ^^^ ^^x|0)(0|xV + ^^^ ^VxlO) . (Olxt-^. (2.15) 

Note that color-octet operators, such as %Ij'^T°'x\^){^\x^T"''^i cire omitted because they can- 
not contribute to matrix elements between color-singlet heavy-quarkonium states. The 
dimension-8 terms that must be added to the lagrangian include 

K™ ™ion),^3 = ^^^^ \ V^^(-tD) ■ ax|0) (0|xt(-f D) . ^^ 



g-EMJ Sq) 1 

M4 2 

^ 9em{'S,) 1 



V'Vx|0)-(0|xV(-iD)2^ + h.c.l + ... . (2.16) 



M4 2 

We have shown only four of the possible dimension-8 terms. In particular, there are terms 
corresponding to each of the operators shown explicitly in ( |2.11| ). The coefficients of the 
operators in ( p.l5| ) and ( p.l6| ) can be computed as perturbation expansions in ^^(M). 



The decay rate of a heavy quarkonium state H into electromagnetic final states (EM) 
can be expressed in a factored form that is analogous to that given in ( f^.l4D for decays into 
light hadrons: 
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r{H ^ EM) = E ^ ^"Ifr-I^"^^ {Hl^^lCUim) (OlxX^(A)li^), (2.17) 

where /C„ and /C^ are products of the unit color matrix, a spin matrix (the unit matrix or 
0"*), and a polynomial in the covariant derivative D and other fields, as in ( p.l5D and ( p.l6[ ). 



The possible electromagnetic final states EM include the multiphoton states 77 and 87 and 
the lepton pairs £~^£~ , where £ = e, fi, r. 

H. Computation of the Coefficients of the 4-Ferniion Operators 

The nonperturbative long-distance dynamics of QCD is described equally well by full 
QCD and by NRQCD. The perturbation expansions for full QCD and NRQCD also give 
equivalent descriptions of the long-distance dynamics, although the description is incorrect. 
For example, perturbation theory allows quarks and antiquarks to appear as asymptotic 
states. However, because the coefficients of the NRQCD operators are insensitive to the long- 
distance dynamics, we can exploit the equivalence of perturbative QCD and perturbative 
NRQCD at long distances as a device to calculate the coefficients of the four-fermion opera- 
tors. We compute in perturbation theory in full QCD the annihilation part A{QQ —>■ QQ) of 
the scattering amplitude for an on-shell quark and antiquark with small relative momenta. 
Then we use perturbation theory in NRQCD to compute the matrix elements of 4-fermion 
operators C„ between on-shell QQ states. The short-distance coefficients are determined by 
the matching condition 



A(QQ -^ QQ) 



= E#^(Wia(A)iQQ) 

pert. QCD n ^^-^ 



(2.18) 

pert. NRQCD 



By expanding the left and right sides of ( |2.18| ) as Taylor series in the relative momenta p 



and p' of the initial and final QQ pairs, we can identify the coefficients of the individual 
operators. These correspond to the infrared-finite parts of the parton-level amplitudes for 
QQ scattering. Because of the equivalence of NRQCD and full QCD at long distances, all 
of the infrared divergences contained in A{QQ -^ QQ) on the left side of ( |2.18| ) reside on 
the right side in the NRQCD matrix elements {QQ\On{^)\QQ) ■ 

30 



The application of the matching condition (|2.18|) is illustrated in Appendix A. The 



imaginary parts of the coefficients /„ that enter into the annihilation rates of S-wave states 
through relative order f ^ and the annihilation rates of P-wave states at leading order in v 
are computed to order a^. In order to illustrate the use of the matching condition ( |2.18| ) 
beyond leading order in a^, we also calculate the coefficient Im/i(^S'o) at next-to- leading 
order in a,. 
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III. MATRIX ELEMENTS FOR HEAVY QUARKONIUM 



The factorization formula ( 2.14|) expresses the decay rate of an arbitrary heavy quarko- 



nium state H into hght hadrons as a sum over all 4-fermion operators On- If we truncate 
the expansion at a given order in the heavy-quark velocity f , then only finitely many of 
the operators contribute. In this section, we show how the number of independent matrix 
elements can be reduced further by exploiting heavy-quark spin symmetry and by using the 
vacuum-saturation approximation. We identify the matrix elements that contribute to the 
decays of S-wave states through relative order v"^ and the matrix elements that contribute 
to the decays of P-wave states at leading order in v. We also discuss the relation between 
these matrix elements and Coulomb-gauge wavefunctions, as well as the dependence of the 
matrix elements on the factorization scale. For the sake of clarity, we use the lowest S-wave 
and P-wave states of charmonium for the purpose of illustration. However, our results apply 
equally well to other sets of S-wave and P-wave states, and they can be extended readily 
to higher orbital-angular-momentum states as well. The lowest-lying S-wave states in the 
charmonium system are the J^^ = 0^^ state rjc and the 1 state J/ip (henceforth referred 
to simply as ip). The lowest-lying P-wave states are the 1^* state he and the J'^^ states 
Xcj, J = 0,1,2. 

A. Powers of Velocity 

We wish to determine the relative importance of the matrix elements {H\On\H) of 4- 
fermion operators On for a heavy quarkonium state \H). The velocity-scaling rules in Table | 
suggest that {H\On\H) is of the same order in v for all the dimension-6 operators in ( |2.9|) , and 
that all the dimension-8 operators in ( p.ll| ) are down by a power of f ^. There can, however. 



be additional suppression by powers of v, depending on the quantum numbers of the state 
H. The velocity-scaling rules in Table I give the correct result only if the operator On 
annihilates and creates a color-singlet QQ pair with the same angular-momentum '^^'^^Lj as 
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the QQ pair in the dominant Fock state \QQ) of the state \H). (In the notation for 4-ferinion 
operators used in (|2.9|) and (|2.11|) , the subscript 1 or 8 and the argument "^^^^Lj indicate 



the color and angular-momentum state of the QQ pair that is annihilated and created by 
the operator.) The matrix element {H\On\H) is suppressed by only one additional power 
of f^, relative to the velocity-scaling rules in Table 1, if On annihilates and creates QQ 
pairs in the same color-spin-orbital state as appears in one of the Fock states \QQg) whose 
probability is of order f ^. In particular, if the dominant QQ component is "^^^^Lj, the QQ 
pair in the component \QQg) must be in a color-octet state with spin quantum number 5* 
and orbital-angular- momentum quantum number L ± 1. The matrix element is suppressed 
by v^ relative to the velocity-scaling rules in Table I if On annihilates and creates QQ pairs 
in the same color-spin-orbital state as appears in one of the Fock states \QQg) that can be 
obtained from the dominant Fock state by a spin-flip transition. In such a Fock state, the 
QQ pair must be in a color-octet state with the same orbital-angular-momentum quantum 
number L as in the dominant \QQ) state, but with different total spin quantum number. 
In all other cases, the matrix element is down by f ^ or more relative to the velocity-scaling 
result from Table I. 

If perturbation theory remained accurate down to the scale Mv, then the spin-flip ma- 
trix elements would be suppressed by an additional power of v. The reason for this is that 
the contribution to a spin-flip matrix element that is suppressed by only v^ relative to the 
velocity-scaling rules is power ultraviolet divergent. Therefore, one could carry out a renor- 
malization of the matrix element in which this contribution is subtracted. The corresponding 
contribution to the decay rate would then reside in the short-distance coefficient of the ma- 
trix element that is associated with the dominant Fock state. (Such a subtraction is carried 
out automatically if dimensional regularization is used to cut off the ultraviolet divergences 
in the matrix element.) Once the subtraction has been made, the leading contribution to 
the spin-flip matrix element comes from the scale Mv"^. It is subject to the usual multipole 
suppression and scales as v'^ relative to the velocity-scaling rules. In practice, one usually 
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makes such subtractions perturbatively. It is not clear, in the charmonium and bottomo- 
nium systems, that perturbation theory is sufficiently accurate at the scale Mv to remove 
the v^ contribution completely. Therefore, we assume in the error estimates below that the 
spin-flip matrix elements scale as v^ relative to the velocity-scaling rules. 

We can now identify the operators that contribute to the annihilation of the rj^ into light 
hadrons through relative order v"^. The J^*" = ^ state \ric) consists predominantly of the 
Fock state \QQ), with the QQ pair in a color-singlet ^5*0 state, but it also contains, with 
probability of order f^, the Fock state \QQg), with the QQ pair in a color-octet ^Pi state. 
The QQ pair in the dominant Fock state \QQ) is annihilated and created by the leading 
operator Oi(^So) = ■?/'^xX^V^? aiid by the operator Vi(^So) = ^p'^xxK"^^)'^^ + ^-c, which 
is down by f^. All other operators are suppressed by v^ or more relative to Oi(^Sq). For 
example, the operator Os(^Pi) = ■?/'^(— |D)T"x-x^(— fD)T"-?/' scales asf^ relative to Oi(^So), 
but it contributes through the Fock state \QQg), which gives an additional suppression by 
v"^. The operator Oi{^Si) = ip'^crx ■ X^^i^ scales as Oi(}So), but it contributes through the 
Fock state \QQgg), and is therefore suppressed by an additional probability factor of f*^. As 
a final example, the operator V(-?/'^x)'^(x^'^) scales as v^ relative to Oi(}So), but its matrix 
element is proportional to the total momentum of the QQ pair, which vanishes for the Fock 
state \QQ) in the meson rest frame. Thus, there are only two operators that contribute to 
the decay rate of the r]c into light hadrons through relative order w^: 

+ ' '"^f ^°^ {Vc\V,CSo)\Vc) + 0{vn (3.1) 

The analysis of the operators that contribute to the decays of the ip is similar to that for 
the r]c- The 1 state {ip) consists predominantly of the Fock state \QQ), with the QQ pair 
in a color-singlet ^5*1 state, but it contains, with a probability of order f^, the Fock state 
\QQg)i with the QQ pair in a color-octet ^Pq) ^-Pi, or ^P2 state. At leading order in f , only 
the operator Oi(^Si) contributes to the decay rate of the ip into light hadrons. At relative 
order f^, the only additional contribution comes from the operator Vii^Si). As in the case 
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of the rjc, all contributions from Fock states containing dynamical gluons, such as \QQg), 
are of order v^T or higher. 

We next determine the operators that contribute to the annihilation decays of the P- 
wave states at leading order in v. In contrast to the S-wave states, we find that Fock 
states containing dynamical gluons play an important role. The l'^ state \hc) consists 
predominantly of the Fock state \QQ), with the QQ pair in a color-singlet ^Pi state, but it 
has a probability of order v^ for the Fock state \QQg), with the QQ pair in a color-octet ^5*0 
or ^D2 state. The Fock state \QQ) is created and annihilated by the dimension-8 operator 
OiCPi) = ^^(-iD)x ■ x^(-|D)V'. The Fock state \QQg), with the QQ pair in a color- 
octet "'^S'o state, also contributes to the decay at the same order in v through the operator 
C^sCSo) = ip''T°'xX^T"''^- The operator scales as f ~^ relative to Oi{^Pi), but there is also 
a suppression factor of f ^ from the probability for the QQg state. Thus, the Fock state 
\QQg), which contains a dynamical gluon, contributes to the decay rate into light hadrons 
at the same order in v as the dominant Fock state \QQ)- The resulting expression for the 
decay rate is 

nhc - LH) = ^ ^"^^['^"^ (/.e|Oi(^Pi)|M 

+ ^ ^"^f ^°^ {K\OsCSo)\K) + 0{v^T). (3.2) 

The analysis of the operators that contribute to the decays of the Xco, Xci, and Xc2 is 
similar to that for he- The J^+ state \Xcj) consists predominantly of the Fock state \QQ), 
with the QQ pair in a color-singlet ^Pj state. It also contains, with a probability of order 
f^, the Fock state \QQg), with the QQ pair in a color-octet ^5*1, ^-Di, ^D2, or ^D^ state. 
The Fock state \QQ) contributes to the annihilation at leading order in v through the 
dimension-8 operator Oi{^Pj). The Fock state \QQg), with the QQ pair in a color-octet ^Si 
state, also contributes to the annihilation rate at the same order in v, through the operator 

The analysis of the electromagnetic operators that contribute to the decays of S-wave 
and P-wave states is identical to that of the operators that contribute to the decays into 
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light hadrons, except that there are no color-octet operators. We find, therefore, that there 
are two operators that contribute to the decay of the rjc into two photons through relative 
order v^: V'^x|0)(0|xV ^^^ V^^x|0)(0|x^(— fD)^'?/' + h.c. Thus, the decay rate for rjc — > 77 
is 

+ ^^""^f^"^ Re {{Vc\^P^x\0){0\x\-lBr^P\Vc)) + 0{v'r). (3.3) 

Here, it is expressed in terms of the vacuum-to-r^c matrix elements (Olx^V^l^c) and 
(0|x^(— fD)^'?/'|^c)- Similarly the decay rate for ip -^ e+e^ can be computed at relative 
order f^ in terms of (0|x^cr?/'|?/') and (0|x^cr(— |D)^^|^). For the decay Xco -^ 77) the 
only operator that contributes at leading order in v is ip'^^—^D ■ cr)x|0)(0|x^(— |D ■ cr)ip, 
so the decay rate can be expressed in terms of the single vacuum-to-Xco niatrix element 
(0|x^(— fD ■ (t)iP\xcq)- Similarly the decay rate for Xc2 -^ 77 can be calculated to leading 
order in v in terms of the matrix element (0|x^(— |-D^V-'^)'?/'|Xc2) only. 

B. Heavy-Quark Spin Symmetry 

Heavy-quark spin symmetry provides approximate relations between matrix elements for 
the various spin states of a given radial and orbital excitation of heavy quarkonium. The 
leading violations of heavy-quark spin symmetry come from the spin-fiip terms in ( p.5D , 
whose effects are of relative order f ^. Consequently, the equalities that follow from heavy- 
quark spin symmetry hold only through relative order v"^. Nevertheless, these relations 
can significantly reduce the number of independent matrix elements that contribute to the 
decays of the various spin states. 

When applied to the S-wave states, heavy-quark spin symmetry relates the rjc and the 
3 spin states of the ip. For the matrix elements that contribute to their decays into light 
hadrons through relative order f ^, the spin-symmetry relations are 

{^lO^CS^m = {v,\O^CSo)\vc) (l + 0{v')], (3.4a) 



36 



{^P\V^CS,m = {vc\V^CSo)\Vc) il + Oiv')j. (3.4b) 

The 12 spin states of the P-wave states h^ Xco, Xci, and Xc2 form a spin-symmetry multiplet. 
The spin-symmetry relations between the matrix elements that contribute to the decays of 
the P-wave states into light hadrons at leading order in v are 

{Xcj\OiCPj)\Xcj) = {K\OiCP,)\K) (i + 0{v^)\ J = 0,1,2, (3.5a) 

{Xcj\OseS,)\xcj) = {K\OsCSo)\K)(l + 0{v')\ J = 0,l,2. (3.5b) 

Heavy-quark spin symmetry also relates the matrix elements that contribute to the 
electromagnetic annihilation decay rates. For the matrix elements that contribute to the 
decays of rjc and ip, the spin-symmetry relations are 

e* • (01x^7^1^(6)) = (OlxVl^c) (l + 0{v^)\ (3.6a) 

e* ■ {0\x^a{-lB)'^\^{e)) = {0\x\-^fiMvc) U + 0{v')\ (3.6b) 

where the polarization vector e satisfies e* ■ e = 1. For the matrix elements that contribute 
to the decays of the Xco and Xc2 into two photons, the spin-symmetry relations are 

e^^-*(0|xH|B(V))^|Xe2(6)) = -^ (Olxn^D-^)^lxeo) (l + 0{v')\ (3.7) 

where the symmetric polarization tensor e*-' satisfies tr(e) = and tr(e^e) = 1. 

C. Vacuum- Saturation Approximation 



The 4-fermion operators in (|2.9| ) and ( p.ll| ) that contribute to the decays of heavy quarko- 



nium into light hadrons are distinct from those in ( |2.15| ) and ( |2.16[ ) that contribute to elec 



tromagnetic annihilation. The electromagnetic matrix elements can be obtained from the 
corresponding light hadronic matrix elements by making use of the "vacuum-saturation ap- 
proximation" : insert a complete set of light hadronic states J2x\-^){-^\ between %"'" and 
X and assume that the sum is "saturated" by the lowest-energy state, the QCD vacuum 
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|0). The vacuum-saturation approximation has been used in many other contexts in parti- 
cle physics, but it is usually just a simplifying assumption without any rigorous basis. In 
the case of heavy quarkonium, we can show that the vacuum-saturation approximation is 
actually a controlled approximation. 

Consider the matrix element of a color-singlet operator of the form 0„ = ■?/'"'" /C^xX^^nV') 
where /C^ and /C,^ are products of a unit color matrix, a spin matrix (the unit matrix or 
cr*), and a polynomial in D and other fields. The vacuum-saturation approximation to the 
matrix element {H\On\H) is obtained by inserting a complete set of states \X) between x 
and x^ 1 and assuming that the sum is well approximated by the term involving the vacuum 
state |0): 

X 

^(if|^t;C;x|0)(0|xX^|i^) (3.8) 

If the last step in ( |3.8| ) is to be a controlled approximation, we must show that the contri- 
butions from all other states, such as multigluon states, are suppressed by powers of v. One 
example for which this can be done is the matrix element {f]c\ip'^ XX^ '^\Vc) ■ I^^ the vacuum- 
saturation approximation, the last hne of ( p3.8D reduces to |(0|xVl'7c)P- This approximation 
would be exact if the rjc were a pure ^Sq QQ state. The point-like operator x^i' would then 
annihilate the r]^ completely, leaving the vacuum state. However, the r/c also has Fock-state 
components, such as \QQg) and \QQgg), which include dynamical gluons. Corrections to 
the vacuum-saturation approximation can be attributed to contributions from intermedi- 
ate states |X) containing such dynamical gluons. For the matrix element (^dV^^XX^V'l^c), 
a single-gluon intermediate state is forbidden by color conservation, so the leading correc- 
tions to the vacuum-saturation approximation come from two-gluon intermediate states \gg)- 
The leading contribution to {ric\'ip^x\gg){g9\x^'^\Vc) comes from the \QQgg) component of 
the rjc, which has a probability of order f ^. Thus, the vacuum-saturation approximation for 
{Vc\'4''^XX^4'\Vc) holds up to corrections of relative order f*^. 

The vacuum-saturation approximation holds up to corrections of relative order v'^ for 
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any matrix element in which the operator creates and annihilates the dominant QQ compo- 
nent of the quarkonium state. For the matrix elements that contribute to the decay rates 
of the S-wave states into light hadrons through relative order f^, the vacuum-saturation 
approximation gives 



{r^,\0,{'So)\Vc) = \{Q\xH\ric) 



1 + Oiv' 



1 + Oiv' 



{ric\Vi{^So)\r]c)= Re ({r]cWx\^){^\xK-h^)MVc)) \l + 0{v 



\V,{'S,)\i,) = Re ((V^I^VxlO) ■ (0|xV(-fD)2^|^)) 1 + 0{v') 



(3.9a) 
(3.9b) 
(3.9c) 
(3.9d) 



In the case of P-wave states, the vacuum-saturation approximation can be applied to the 
matrix elements of the color-singlet 4-fermion operators of dimension 8: 



(/ie|Oi(^Pi)IM= (oix^(-fD)^lM 



1 + Oiv' 



(Xeo|0i('Po)|Xco) = - (0|xH-fD ■ <T)ij\Xc0) 



1 + 0(v' 



{Xci\0,{'Pi)\Xci) = - |(0|xH-fD X a)^\xci) 1 + 0{v')], 



{Xc2\0,eP2)\Xc2) = E kO|X^(-tS(V^-))^|Xe2) 



1 + 0(v' 



(3.10a) 
(3.10b) 
(3.10c) 
(3.10d) 



v 



The vacuum-saturation approximation cannot be applied to matrix elements of color-octet 
operators, such as {hc\ip'^T°'xx'^T°'ilj\hc) , because the matrix element {X\x''T""ip\hc) vanishes 
if (X| is the vacuum or any other color-singlet state. 

D. Relation to Wavefunctions 

In most previous work on the annihilation decays of heavy quarkonium, the nonperturba- 
tive factors in the decay rates were expressed in terms of wavefunctions, or their derivatives, 
evaluated at the origin. These "wavefunctions" were often identified with the Schrodinger 
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wavefunctions calculated in potential models for heavy quarkonium. The wavefunction fac- 
tors were never given rigorous field-theoretic definitions, so the accuracy of the approxima- 
tions that were involved was always vague. By expressing the decay rates in terms of matrix 
elements of NRQCD, we have provided a rigorous field-theoretic definition of the nonper- 
turbative factors in the decay rates. Since heavy quarks and antiquarks are described in 
NRQCD by a Schrodinger field theory, the nonrelativistic wavefunctions can also be given 
rigorous field-theoretic definitions, and their relations to the nonperturbative factors in the 
decay rates can be clarified. 

Nonrelativistic Coulomb-gauge wavefunctions can be defined naturally as NRQCD 
Bethe-Salpeter QQ wavefunctions, evaluated at equal time. For example, the radial wave- 
function Rr^^{r) for the rjc can be defined as 



V47r v^^Vc 



(3.11) 

Coulomb 

The Pauli spinor fields ■?/'(r/2) and %''"(— r/2) are understood to be evaluated at the same 
time t = 0. The factor l/^/ivr on the left is the spherical harmonic Yoq^t), while the factor of 
y/2Nc on the left takes into account the traces of the spin wavefunction 5™-+™ '°/v^ and the 
color wavefunction 5"^^ / ^/Wl of the \QQ) component of the rj,.. In the absence of a regulator, 
the wave function or its derivatives may be singular as r ^ 0. We can define regularized 
"radial wavefunctions at the origin" R^^{A) for rjc and R^{A) for ip by 

R^^(A)^^{0\x^^(A)\r^^), (3.12a) 

V -' 'c 

2^ 




R4A) e ^ ^/^ (OlxVV^(A)lV^(e)) , (3.12b) 

where e is the polarization vector of the ip. The local operators x^V'(^) ^-iid x^f^'^Pi^) can be 
defined by dimensional regularization with scale A, together with minimal subtraction. They 
can also be defined by a lattice regulator, or any other convenient regularization scheme. As 
is suggested by the overline, the intuitive interpretation of i?^^(A) and -R^(A) is that they 
are the radial wavefunctions averaged over a region of size 1/A centered at the origin. Note 
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that the regularized matrix elements in ( |3.12|) are precisely those that enter into the decay 



rates for for rj^ —>■ 77 and tp —>■ e+e~ at leading order in v. 

The matrix element {0\x^ {— ^J^Yiplrjc) that contributes to the decay rate for rjc —>■ 77 at 
relative order f^ can also be related to the Coulomb-gauge wavefunction defined in ( p. 111 ). 
By the velocity-scaling rules of Table |I|, it differs from the matrix element (0|x^(— f V)^'?/'|r7c) 
in Coulomb gauge only at relative order f^. With appropriate regularization, the latter 
matrix element can be identified with the limit as r — i> of —V R{r), where R{r) is the 
radial wavefunction defined in (|3.11|) . The operator x^ll^)^"^ contains a linear ultraviolet 



divergence proportional to xV(A), which we subtract, and a logarithmic divergence that is 
cut off at the scale A. This subtraction and cutoff define a renormalized laplacian of the 



radial wavefunction at the origin, which we denote by V^-R^^: 

(3.13) 



V^i?,.(A) = \h7- {0\xHl^r^{^)\Vc) 



Coulomb 



The analogous quantity V^R^ for the ip can be defined in a similar way. The corresponding 
gauge-invariant matrix elements differ from V^-R^^ and V^i?^ only at relative order f^: 

(0|xn-fD)2^(A)|r/,) = -J^ V^(A) U + 0{v')\. (3.14a) 




xV(-fD)2^(A)|^(e)) = - \-^ V^R4A) 1 + 0{v')]. (3.14b) 



The intuitive interpretations of V^-R,,^(A) and V^-R^(A) are somewhat obscured by the 
subtractions needed to define the renormalized matrix elements. 

Heavy-quark spin symmetry implies that the wavefunctions of the rjc and ip are identical 
up to corrections of relative order v"^: 

R^r) = R^Xr) fl + 0{v')\. (3.15) 

It is convenient to introduce an average radial wavefunction Rs{r) for the IS states 1]^ and 
ip, which can be used when the differences of relative order f ^ can be neglected: 
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Because of the heavy-quark spin symmetry, the regularized quantities i?^^(A) and -R^(A) 



differ only at relative order f^, as do the renormalized quantities V^-R,j^(A) and V^-R^(A). 



Weighted averages Rs{-^) and V^Rs{A) for the S-wave states can be defined as in ( p.lGI) . 
The S-wave radial wavefunction that is computed in nonrelativistic potential models can 
be interpreted as a phenomenological estimate of the wavefunction ( p.l6| ). Thus, the value 
Rsir = 0) that is obtained from potential models can be used as an estimate of the reg- 



ularized quantity Rs{A.) at a scale A of order Mv. The relation between V^Rs{A.) and 
'V'^Rs{r = 0) in potential models is more obscure, because of the subtraction that is re- 
quired to define the renormalized matrix element in ( p. 131) , and because V Rs{r) diverges 



linearly as r ^ if the potential is Coulombic at short distances. 

Nonrelativistic wavefunctions for the P-wave states can be defined through matrix ele- 



ments in Coulomb gauge that are analogous to (|3.11| ). For example, the radial wavefunction 
Rhc.ij') for the he can be defined as 



Rj^Xr) h/|^r-el ^ -^ (0|x^(-r/2) ^(+r/2)|M6)) 



(3.17) 

Coulomb 



where the polarization vector satisfies e ■ e* = 1. A regularized derivative of the radial 



wavefunction at the origin i?^^(A) can be defined by 

(3.18) 



i?',^(A)6 ^ J^ (0|xn|V)^(A)|/^c(e)) 



Coulomb 



Analogous quantities R' (A) can be defined for the XcJ states. Since no subtractions are 



required in order to define the operator on the right side of (|3.18|) , the quantity i?^^(A) 



has a straightforward intuitive interpretation as the derivative of the radial wavefunction 
averaged over a region of size 1/A centered at the origin. In the gauge-invariant analog of 
the matrix element ( p.l8|) , the derivative V is replaced by the covariant derivative D. From 



the velocity-scaling rules of Table |I|, we see that these matrix elements differ only at relative 
order v^: 



(0|xniD)^(A)|/.,(6)) = f-^ RliA) e (l + 0{v')\ (3.19a) 
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i= (Olxn^D ■ ^)^(A)IXco) = \/^ R',JA) (l + 0{v')\ (3.19b) 



i= (Olxn-fD X ^)^(A)|xei(6)) = ^^ i?;^^(A) 6 (^1 + Oiv^yj, (3.19c) 

(0|xn|S(V^))^(A)|Xc2(e)) = ^^^A^) ^'' (l + 0(t;^)). (3.19d) 



By heavy-quark spin symmetry, i?^^(A) differs from i?' ^(A), J = 0,1,2, only at relative 
order v^ . For applications in which v^ corrections can be neglected, these wavefunctions can 
all be replaced by the average over the 16 P-wave spin states, which we denote by -Rp(A). 
The value -Rp(O) for the derivative of the radial wavefunction at the origin that is obtained 
from nonrelativistic potential models can be interpreted as a phenomenological estimate of 
the regularized quantity -Rp(A) at a scale A of order Mv. 

The vacuum-saturation approximation discussed in Section |111C| allows the matrix ele- 
ments of some 4-fermion operators to be expressed in terms of the regularized and renor- 
malized wavefunction parameters defined above. Combining (|3.9|) with ( 3.12 ) and ( p. 141) , 



we obtain the following expressions for the matrix elements that contribute to the decays of 
the r]c and the ?/> into light hadrons: 

(r/c|Oi('5o)|r/,) = ^\^^{\^ 0{v^)\. (3.20a) 

(^|Oi(35i)|^) = §N' (l + 0{y')\ (3-20b) 

(r7c|Pi('5o)|r/,) = -^ Re(:R^* V^) [\ + O^v")]^, (3.20c) 



(^|Pi(^5i)|^) = -^ Re(i?/ V2/?5) 1 + 0{v')\. (3.20d) 



In ( |3.20c| ) and (|3.20d|) , we have used heavy-quark spin symmetry to replace i?^^ and K 



"i/j 



by their weighted average Rs and to replace V^-R.^^ and V^-R^ by V^-Rs without any loss 
of accuracy. If we were to make the same replacement in ( p.20a| ) and (|3.20b|) , the relative 
accuracy would be decreased to v^ . For the decays of the P-wave states into light hadrons 
at leading order in f , the vacuum-saturation approximation together with heavy-quark spin 
symmetry can be used to express all the color-singlet matrix elements in terms of the average 
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regularized quantity R'p. Combining ( |3.10| ) with ( p.l9| ), we obtain the approximations 



{h,\0,{^P,)\K) = ^|i?pf (l + 0{v^)\ (3.21a) 



3iV, 



2 



{x,j\0^{:'Pj)\Xcj) = ^R'p [I + 0{v')l J = 0,l,2. (3.21b) 



E. Factorization-Scale Dependence 



The matrix elements {H\On\H) that appear in the factorization formula ( |2.14| ) are ultra- 
violet finite only if the 4-fermion operators On are properly regularized. The regularization 
introduces dependence on the ultraviolet cutoff A of NRQCD, and this cutoff-dependence 
must be understood in order to make quantitative predictions. We assume that the operator 
On is normal-ordered: (0|(9„|0) = 0. This guarantees that in the matrix element {H\On\H), 
the operator On annihilates the heavy quark and antiquark in the initial quarkonium state 
\H). In addition to normal-ordering, regularization is needed to control power and logarith- 
mic divergences. If a cutoff A is imposed on loop momenta, there are power divergences in 
{H\On\H) that are proportional to A^, where p = 1, 2, . . .. If the operator On has dimension 
dni then the coefficient of A^ is, by dimensional analysis, the sum of matrix elements of 
4-fermion operators of dimension dn—p or larger. If the dimension is larger than dn — p, the 
extra dimensions are balanced by powers of 1/M. Similarly the coefficients of logarithmic 
divergences are proportional to matrix elements of 4-fermion operators of dimension dn or 
larger. 

The power and logarithmic divergences associated with loop corrections to NRQCD 
operators can be regularized by a variety of means. A convenient regularization scheme 
for analytic calculations is dimensional regularization with minimal subtraction. The scale 
associated with the dimensional regularization then plays the role of the NRQCD cutoff A. 
An advantage in using a mass-independent regulator, such as dimensional regularization, is 
that power divergences are automatically discarded. In other approaches, such as lattice 
regularization, the regularized operator may contain divergences that are proportional to 
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powers of A. These power divergences are simply artifacts of the regularization scheme and 
have no physical content. Since physical quantities are renormalization-group invariants, 
they have no dependence on A. Hence, any power divergences in NRQCD operator matrix 
elements must ultimately be cancelled by power divergences in operator coefficients. 

Once one has removed the power divergences, either by employing a mass-independent 
regularization scheme or by making explicit subtractions, the 4-fermion operators satisfy 
simple evolution equations of the form 

^^0„W = J:^0,(A). (3.22) 

k 

where the sum ranges over all 4-fermion operators Ok{A) with dimensions dk > (i„. The 
anomalous-dimension coefficients 7„fc(A) are computable as power series in the running cou- 
pling constant ^^(A). For dn = dk, the coefficients 7„fc are of order a^, because logarithmic 
ultraviolet divergences in one-loop diagrams in NRQCD arise only from transverse gluons, 
whose coupling to the heavy quark lines brings in a factor of v. The coefficients jnk for 
dn = Q and dk = 8 are computed to order as in Appendix B. 

By taking the matrix elements of ( |3.22| ) between heavy quarkonium states \H), we obtain 
the evolution equations for the matrix elements {H\On{A.)\H) that appear in the general 
factorization formula (|2.14| ): 



AUnmAm = j:^^iH\OkiA)\H). (3.23) 

k 

The leading v behavior of the matrix elements can be determined by using the velocity- 
scaling rules developed in the previous sections. At any given order in v, there is only 
a finite number of terms that contribute to the evolution equation ( p.23|) . The evolution 



equations for the dimension-6 4-fermion operators are calculated to order a^ in Appendix B. 
The operator evolution equations for Oi(^So) and Oi{^Si) are given in ( [B16| ) and ( [BlQal ). 



Taking the matrix elements of these equations and keeping only those terms on the right 
sides that are of relative order f ^, we find that only the operators Vi survive, and we obtain 
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A^(r7c|0i(^5o)|r^.) = -^§^^(^c|Pi(^5o)|r/e), (3.24a) 

K^{i,\0,{'S,m = -^-^^i^l-PieS.M), (3.24b) 

where Cp = (N^ — l)/{2Nc). If the evolution equations are truncated at leading order 
in V, the right sides of ( ^.24aD and ( pj.24b| ) vanish and we find that the matrix elements 
(?7c|Oi("'^S'o)|77c) and {il!\Oi{^Si)\il!) are renormalization-scale invariant through order a^. The 
dimension-8 matrix elements {i]c\'Pi(^So)\t]c) and {ip\Vi{^Si)\il!) are also renormalization- 
scale invariant through order a^ and at leading order in v. 

^^{Vc\riCSo)\Vc) = 0, (3.25a) 

A^{^\V,eS^m = 0. (3.25b) 

Truncated at order a^, the evolution equations can be solved analytically for the A- 
dependence of the matrix elements. For example, the solution to ( p.24a| ) is 

{Vc\O^CS,;A)\vc) = {vMCSo;Ao)\Vc) - ^^ l^S (^^) {Vc\V,CSo)\Vc) , (3.26) 

where /3o = {llNc — 2nf)/Q is the first coefficient in the beta function for QCD with Uf 
flavors of light quarks: fi{d/dfi)as = — /^oa^/vr + . . .. 

We next consider the evolution of the matrix elements that contribute to P-wave 
annihilation at leading order in v. The color-singlet dimension-8 matrix elements are 
renormalization-scale- invariant to this order in a^: 

A^{h,\OiCPi)\hc)=0, (3.27a) 

A^{Xcj\OiCPj)\xcj)=0, J = 0,1, 2. (3.27b) 

Taking the matrix elements of ( |B17] ) and ( [B19b| ) in Appendix B, we find that the color-octet 
dimension-6 matrix elements have nontrivial scaling behavior at order as'- 

A^{K\OsCSo)\K) = ^]^(^e|Oi(^Pi)IM, (3.28a) 

A^{Xcj\OseS,)\xcj) = f^2^^2 {Xcj\0,CPj)\xcj). (3.28b) 
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To this order in a^, we find that the evolution equations can be solved analytically. For 
example, the solution to (|3.28a|) is 



{K\OsCSo;A)\h^) = {K\OsCSo;Ao)\h.) + J^^^ log (^y) {h^\0,{' P,)\K) . 



(3.29) 



The solution ( p.29|) to the evolution equation can be used to provide a crude estimate of 



the color-octet matrix element (/idOsl^'S'o; A)|/ic) in terms of the color-singlet matrix element 
{hc\Oi{^ Pi)\h^ . Suppose that we approximate ( p.29|) by the evolution term on the right side. 
The evolution term is largest, relative to the matrix element (/idOsl^'S'o; Ao)|/ic), when the 
scales Aq and A are as widely separated as possible. However, the logarithmic evolution 
holds only down to scales of order Mv. Thus, we choose Aq = Mv. Then, setting A = M, 
neglecting the initial matrix element in ( p.29| ), and assuming that ^^(Ao) = as{Mv) ~ f , 



we find that ( p.29| ) reduces to 



ACf 



i^'\0'('s,.M)\K) « 5^^^ log [^^^) {Km^mK) . (3.30) 

The same method can be used to obtain crude estimates for the corresponding matrix 
elements for the XcJ states: 

The terms that we have retained in obtaining these estimates are enhanced by one power 
of log[f/Q;s(M)] relative to the terms that we have neglected. Since this is not a large 
enhancement factor, particularly in the case of charmonium, these estimates should be 
regarded as giving only the orders of magnitude of the matrix elements. 
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IV. ANNIHILATION DECAYS OF HEAVY QUARKONIUM 



In Section |Tl|, we derived a factorization formula (|2.14| ) for the decay rates of heavy 



quarkonium states into hght hadrons. In this section, we apply this formula to the decays 
of S-wave states, up to corrections of relative order v^, and to the decays of P-wave states, 
up to corrections of relative order v"^. We also treat the decays into the electromagnetic 
final states by using the analogous formula ( |2.17|) . As in Section |n|, we use the lowest-lying 



S-wave and P-wave states of charmonium for the purpose of illustration. 

A. S-wave Annihilation 

Most previous treatments of the annihilation rates of the S-wave states of heavy quarko- 
nium have been restricted to leading order in v. In these analyses, long-distance effects were 
absorbed into a nonperturbative factor [/^^(O)^, where Rs{0) is the radial wavefunction at 
the origin. We improve on these previous treatments by providing a rigorous definition of the 
nonperturbative factor in terms of matrix elements of NRQCD. We also extend the analysis 
of the decay rates to relative order f ^, and show that 3 independent nonperturbative factors 
are sufficient to calculate all the S-wave annihilation rates through this order. 

We first consider the decays of the J^*" = ^ state r/c and the 1 state tp into light 



hadrons. As was shown in Section |III A| , there are only two operators that contribute to 
each of these decay rates through relative order v^. According to ( p.l4|) , the decay rates into 
light hadrons are therefore 

nnc -> LH) = ^^^^^^ iVclO^CSoU) + ^^^^fl^ {Vc\V,CSo)\v.) 

+ 0{v^T), (4.1a) 

r(^ - LH) = ^^^^^^ {^wsm + ^^^^fp^ {^pw^csm 

+ 0{v^r). (4.1b) 

The imaginary parts of the coefficients in ( |4.1| ) are calculated at order a^ in Appendix |A 2| , 
and Im /i(^S'o) and Im fi{^Si) are given through next-to-leading order in a^ in Appendix |A 3 
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According to the factorization formula ( p.l7| ) for electromagnetic annihilation, the decay 
rates for tjc -^ 77 and ip — > e"^e~ are 

2Im/^^(i5o) 



r(r/c -> 77) 
+ 



M2 
2 Im5f^^(iS'o) 



(oixVl^c) 



M4 



Re ((r/e|^t^|o)(o|x^(-|D)V|r7c)) + 0{v^V), 



(4.2a) 



e e 



2Im/ee(=^5i: 



M2 



(OlxV^lV') 



+ 



2Im^ee('>Si 



Re ( (VlV^VxlO) ■ (0|xV(-|D)V 



0(t;^r). (4.2b) 



The decay rate for ip -^ 777 is given by an expression that is identical to ( |4.2b| ), but 
with coefficients fz-fi^Si) and (y'3^(^S'i). The imaginary parts of the coefficients in ( ^.2| ) are 
calculated at order a^ in Appendix |A4| , and order-a^ corrections are given for Im/^^(^5'o), 
Im/ee(^>S'i), and Im/3^('^5'i). The matrix elements in ( [4.1[ ) and ( |4.2| ) can be computed using 
lattice simulations of NRQCD. Since matrix elements of relative order v^ have been omitted, 
there is nothing to be gained by computing the dimension-6 matrix elements to an accuracy 
of better than v"^. Similarly, the dimension-8 matrix elements need be computed only at 
leading order in v. 

At the level of accuracy in ([4.1|) and ( [4.2|) , the matrix elements are not all independent. 
The vacuum-saturation approximation (|3.9| ) can be used to express the 4-fermion matrix ele- 
ments in ( [4.1| ) in terms of the vacuum-to-quarkonium matrix elements in ( [4.21 ). Furthermore, 
the heavy-quark spin-symmetry relation ( p.6| ) can be used to equate the matrix elements in 
the second terms on the right sides of ( |4.2a| ) and ( |4.2b| ). The net result is that the 8 matrix 
elements in ( |4.1| ) and ( [4.2| ) can be reduced to 3 independent nonperturbative quantities. 



which we can take to be |-R.^J^, l-R^'P' ^"^^ Re(i?5* V^Rs)- The resulting expressions for the 
decay rates are 



T{7ic ^ LH) 
r(z^ -^ LH) 

r(r?c -^ 77) 



7rM2 

7rM2 
A^cIm/^T,(^S'o; 



R 



Vc 



2 N.Img^CSo) ^ ,—, 



TlM^ 



Re{Rs*Vms) + 0{v'T) 



R11 



2 AT, Im (71 (3^1) ^ ,_ 



TXM^ 



Re{Rs*Vms) + 0{v'T) 



(4.3a) 
(4.3b) 



TC 



M2 



Rr, 



2 N^lmg^^i^So) ^-^ 



ttM^ 



Re{Rs*Vms) + 0{v^T), (4.3c) 
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T{il) -^ e+e') 



N.lmfU'S, 



R 



2 N.lmgU^Si) ^ ,^, 



■kM^ 



Re{Rs*Vms) + 0(wT). (4.3d) 



7rM2 
The quantities i?^^, R^, and V^Rs are defined in Section |IIID| in terms of vacuum-to- 



quarkonium matrix elements in NRQCD, and can, tlierefore, be calculated using nonpertur- 
bative methods, such as lattice simulations. They can also be estimated using the wavefunc- 
tions that are obtained from nonrelativistic potential models of quarkonium. Alternatively, 
since there are more decay rates than there are parameters, they can be treated as purely 
phenomenological parameters, to be determined by experiment. 

The approximations to the matrix elements that have been made in (|4.3|) imply restric- 
tions on the order in as{M) to which the coefficients can be included meaningfully. Because 
of the identification of v with as{Mv) in ( p.l|) , we should consider v to be less than or of 
order as{M). There is no point in calculating the coefficients to relative order a" unless we 
have included all operators whose matrix elements are of relative order v"' or less. Hence, 
there is no gain in accuracy if the coefficients of |-Rr;cP ^^^ l-^vP ^^^ calculated beyond 



relative order a^, or if the coefficients of Re(i?5 V^Rs) are calculated beyond relative order 



a... 



If we require accuracy only to leading order in v, then the decay rates in (^4.3|) can be 
simplified further. The difference between i?,,^ and R^ is of relative order v"^, so both can be 



replaced by their weighted average Rs- The factor Re(i?5* V^-R^) is of order f ^ relative to 
|i?5p and can therefore be set to 0. We thereby recover the familiar factorization formulas 
assumed in previous work: 



r{vc -^ LH) 
r(^ -^ LH) 

r(r/c -^ 77) 



iVelm/iCgo) 
7rM2 

7rM2 

7rM2 

7rM2 



R, 



R, 



R, 



R, 



+ 0{v'T), 



(4.4a) 
(4.4b) 
(4.4c) 
(4.4d) 



Because corrections of relative order f2 have been neglected in ( ^.41 ), there is no point in 
calculating the regularized wavefunction at the origin Rs to an accuracy of relative order v"^. 
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Similarly, because of the identification of v with as{Mv)^ there is no increase in accuracy 
if the coefficients of \Rs\^ in ( [4.4| ) are calculated beyond next-to-leading order in as- The 



effects of matrix elements of relative order v"^ are probably more important than perturbative 
corrections to the coefficients that are of relative order a^. 

B. P-^vave Annihilation 

In most previous work on the annihilation decays of P-wave states, it was assumed that 
long-distance effects could be factored into a single nonperturbative quantity |i?p(0)p, where 
i?p(0) is the derivative of the radial wavefunction at the origin. By explicit calculation, 
Barbieri et al. [^H found that the coefficients of |i?p(0)p depend logarithmically on an 
infrared cutoff on the energies of the final-state gluons. In subsequent phenomenological 
applications of these calculations, the infrared cutoff has been identified with the binding 
energy of the quarkonium state, which is of order Mf ^, the inverse of the radius of the bound 
state, which is of order Mf , or the inverse of the confinement radius, which is of order Aqcd- 
It should be clear, however, that the infrared divergence is a signal of the breakdown of the 
factorization assumption upon which the calculation is based. The solution to the problem of 
infrared divergences in the calculation of the P-wave decay rates into light hadrons was first 
presented in Ref. 0]. We will review the resolution of this problem later in this subsection. 



As was shown in Section [III A| , there are two 4-fermion operators that contribute to the 
decay rates of any of the P-wave states into light hadrons at leading order in v. According to 
our factorization formula ( p.l4|) , the decay rates of the four P-wave states into light hadrons 



are 



nhc - LH) = ^ ^"^['^'^ {K\O^CP^)\K) 



+ ^^""^f^"^ {K\OsCSo)\K) + 0{v'T), (4.5a) 



nxcj - LH) = ^ ^"^r^^^ {xcj\o,epj)\xcj) 



1^^ {x,J\OseS^)\xcJ) + Oiv'T), J = 0,1,2. (4.5b) 
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The imaginary parts of the coefficients /i(^Po)7 fi{^P2), fsi^So), and fs{^Si) are calculated 
in order a^ in Appendix [A2|. The color-octet matrix elements in the factorization formulas 



]5|) represent contributions to the annihilation rates from the Fock states \QQg)- Thus, 
we see that the decays of P-wave (and higher orbital-angular-momentum states) can probe 
components of the meson wavefunction that involve dynamical gluons. For the decays of the 
Xco and Xc2 into two photons, there is only one operator that contributes at leading order 



m v: 



2Im/^^(3Po) 1 I,-. ., .- , . ,2 



3 , _ . . . 

2t 



r(Xco - 77) = '-^t^ 3 |(0|xn-fD ■ ct)^\xco) + Oiv^T), (4.6a) 



'TT^ "^ ± \/n\.Ai 
r(Xc2 - 77) = ^^"^tf^'^ E |(0|xn-tS(V^))^|Xe2) ' + Oiv'T). (4.6b) 






The coefficients Ira /^^(^Pq) and Im/^^(^P2) are calculated at order a^ in Appendix |A4 



and the order-a^ corrections are given as well. There is no increase in accuracy if the matrix 
elements in ( [4.5|) and (|4.6| ) are calculated to an accuracy of relative order v"^, since matrix 



elements of relative order v have been omitted. Because of the identification (|2.1|) of v with 



as{Mv), there is no increase in accuracy if the coefficients in ( [4. 5] ) and ( [4. 61 ) are calculated 



beyond next-to- leading order in as{M). Perturbative corrections of relative order a^(M) 
are probably less important than contributions of other matrix elements of relative order v"^. 
To the order in v that is being considered in ( [4.5|) and (|4.6| ), the matrix elements are 
not all independent. The vacuum-saturation approximation ( p.lO|) can be used to express 
the matrix elements of Oi{^Pi) and Oii^Pj) in ([4.5|) in terms of vacuum-to-quarkonium 



matrix elements. These matrix elements can be related to regularized derivatives of radial 
wavefunctions at the origin by using (|3.19|) . Because of the heavy-quark spin symmetry. 



the derivatives of the radial wavefunctions at the origin can all be replaced by the average 
value R'p for the 12 spin states of the P-wave system, without any loss of accuracy. The 
heavy-quark spin-symmetry relation ( |3.5b| ) also implies that the matrix elements of OsCSq) 
and Osi^Si) in ( |4.5| ) are the same, up to corrections of relative order v"^. Thus, the decay 
rates ( |4.5| ) and ( [4. 6] ) can all be expressed in terms of two nonperturbative quantities |-Rpp 
and {hc\Os(^SQ)\hc) (or, alternatively, the average of the color-octet matrix elements for the 
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12 P-wave spin states): 



7rM4 



M2 



r(x. ^ LH) ^ ^-^ ■!■ {f ^) fif ^ ^ iiHiM^i) (.J0.C5„)|M 



ttM^ 



M2 
+ 0(i;2r), J = 0,1, 2, (4.7b) 



3iV,Im^^(^Pj; 



2 



r(XcJ - 77) = "^'^7^7^^^^ ^P ' + 0{v'T), J = 0, 2. (4.7c) 

Since R'p is proportional to a vacuum-to-quarkonium matrix element, it can be calculated 
more easily in lattice NRQCD simulations than can {hc\Os(^ SQ)\hc) , which is a matrix 
element between quarkonium states. 

As we have already mentioned, in the calculations of Barbieri et al. of the P-wave decay 
rates into light hadrons 0,^,^ , a logarithmic dependence on an infrared cutoff appeared in 



the coefficients of |-Rp(0)p. We now explain why this infrared-cutoff dependence is absent 
in the factorization formulas ( [4.7| ). The coefficients of \R'p\^ in ( |4.7| ) depend logarithmically 
on the NRQCD cutoff A. In these coefficients, A plays the same role as did the infrared 
cutoff in the Barbieri et al. calculations. According to the evolution equation ( |3.28aD , 



the matrix element {hc\Os{} SQ)\hc) , also depends logarithmically on A. In this case, A 
plays the role of an ultraviolet cutoff. Because physical quantities, such as decay rates, 
are renormalization-group invariants, the A-dependence in {h(.\0^(}SQ)\h^ cancels the A- 
dependence in the coefficients of |-Rpp in (|4.7| ). Thus, we see that the inclusion of the 
color-octet term proportional to (/ic|08(^'S'o)|/ic) in the factorization formulas removes the 
dependence of the decay rate on an arbitrary infrared cutoff. 

The factorization formulas ( ^.7| ) for the annihilation decays of P-waves at leading order 
in V were first given in Ref. |]10| in the form 



T{h^ -^ LH) = Hi fi{QQ{^Pi) -^ partons) + H^ f^{QQ{^SQ) -^ partons), (4.8a) 

T{xcj -^ LH) = Hi TiiQQCPj) ^ partons) + Hs TsiQQCSi) -^ partons), 

J = 0,1, 2, (4.8b) 

r(Xcj ^ 77) = Hi TiiQQi'Pj) ^ 77), J = 0, 2. (4.8c) 
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The coefficients Fi and Fg in ( [4 .81) are proportional to the annihilation rates of on-shell QQ 
pairs in color-singlet P-wave and color-octet S-wave states, respectively. The nonperturba- 
tive parameters Hi and H^ that were introduced in Ref. ||10[ can be defined rigorously in 
terms of matrix elements in NRQCD: 

Hi = ^^{K\OiCPi)\K), (4.9a) 

Hs{A) = ^{K\OsCSo)\K). (4.9b) 



The factors of 1/M^ and l/M^ in (|Oal ) and {]^^ were chosen in Ref. ^ so that Hi 



and Hs would be the combinations of NRQCD matrix elements and quark masses that are 
determined in experimental measurements of the P-wave decay rates. 



In retrospect, the choice made in Ref. ||T0[ to include factors of 1/M in the definitions of 



Hi and Hg in ( [4.9|) was unfortunate. The factors of 1/M are more properly associated with 
the coefficients Fi and Fg, since they involve short-distance physics at scales of order 1/A 
or less. The factorization formulas ([4.7|) are, therefore, to be preferred over the forms ( [4.8D , 
because they incorporate all effects of the short distance scale 1/M into the coefficients, 
leaving matrix elements that depend only on physics at length scales l/{Mv) and longer. 
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V. PERTURBATIVE FACTORIZATION 

In this section, we sketch the connection between the NRQCD approach and conven- 
tional perturbative methods for demonstrating the factorization of cross sections involving 
large momentum transfer in QCD. In perturbative proofs of factorization, the aim is to 
demonstrate that, to all orders in perturbation theory, infrared and coUinear divergences 
either cancel or can be absorbed into well-defined nonperturbative long-distance quantities. 
Some familiar examples of such nonperturbative quantities are parton distributions in the 
case of deep-inelastic lepton-hadron scattering and fragmentation functions in the case of 
inclusive hadron production at large transverse momentum in e~^e~ annihilation. The cross 
sections can be written as sums of products of long-distance quantities with infrared-safe 
{i.e., short-distance) parton-level cross sections. Our factorization formula for heavy quarko- 
nium annihilation rates is also of this form, and it is illuminating to see how it could be 
derived from a more conventional perturbative analysis. 

A. Topological Factorization 

We remind the reader that, in QCD, infrared (or soft) divergences are logarithmic and 
arise only from the emission of a gluon for which all components of the 4- momentum are 
small. Collinear divergences (or mass singularities) are also logarithmic, and arise when one 
parton (gluon or light quark) splits into two or more partons and all of their 4-momenta are 
parallel. Collinear divergences are cut off by quark masses, which necessarily introduce a 
non-parallel component into the 4-momenta of the splitting partons. 

Let us focus first on the infrared divergences that arise in the radiation of gluons from 
final-state partons and on the collinear divergences that arise in the splitting of a final-state 
parton into collinear partons. The Kinoshita-Lee-Nauenberg theorem [|1^ guarantees that 



all such divergences cancel when one sums over those final-state cuts of a given diagram 
that contribute to an inclusive cross section. For example, the diagram shown in Fig. |l| 
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has three cuts that correspond to gluonic final states, and each cut contains infrared and 
colhnear divergences. However, these divergences cancel when one adds the real-emission 
cut of Fig. |I](b) to the virtual-emission cuts of Figs, ^(a) and |l|(c). This step in the proof 
of perturbative factorization is related to the localization of the annihilation vertex, which 
was discussed in Section lllE, 



Next, let us consider the radiation of gluons from the heavy-quark lines. Such contribu- 
tions are protected from collinear divergences by the heavy-quark mass, so we need consider 
only the possibility of infrared divergences. One key to analyzing the infrared divergences is 
the concept of a "controlling momentum". The essential idea is that the infrared divergence 
associated with an integration over propagators and vertices in some portion of a Feynman 
diagram is cut off by the largest external momentum that enters the propagators. For exam- 
ple, an infrared divergence could potentially arise from the square of the diagram in Fig. ^(c) 
when all components of the 4-momentum of the middle gluon become small. However, be- 
cause of simple kinematics, the other two final-state gluons must both carry large momenta, 
some of whose components are of order M . That large momentum must fiow through the 
heavy-quark propagator to which the soft gluon attaches, and, consequently, it cuts off the 
potential infrared divergence. 

In this example, and in general, the concept of a controlling momentum tells us that an 
infrared divergence can never arise from a soft gluon that attaches to a propagator that is 
off-shell by order M. That means that the infrared- divergent part of a Feynman diagram 
can always be separated from the "short-distance part" by cutting through heavy-quark 
propagators that are off the mass shell by amounts that are much less than M. (By the 
short- distance part, we mean that part of the diagram that includes the hard final-state 
partons and all propagators that are off-shell by order M.) This "topological factorization" 
is the crucial step in a perturbative demonstration of factorization. It implies that the 
infrared divergences can be disentangled from the short-distance part of the diagram and 
absorbed into the long-distance part of the diagram, which also includes the quarkonium 
wavef unctions. 
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The topological factorization of the annihilation rate of heavy quarkonium is represented 
schematically in Fig. 0. The shaded ovals represent the wavefunction for a quarkonium state. 
A typical Fock state contains a QQ pair and zero or more gluons or light quark pairs. The 
short- distance part of the annihilation rate is represented by the circle labelled H (for hard). 
At leading order in w , the only lines that attach to the short-distance part are the incoming 
Q and Q and the outgoing Q and Q. The long-distance part includes the wavefunction of 
the initial meson and its complex conjugate. These wavefunctions are connected by any 
extra partons that may be present in the Fock state, which are represented in Fig. ^ by 
gluon lines. The long-distance part also includes soft-gluon interactions between the extra 
partons, which are represented by the circle labelled S (for soft). 

Once topological factorization has been demonstrated, two additional steps are required 
in order to complete the proof of perturbative factorization. First, one must decouple the 
relative 4-momentum p of the heavy quark and antiquark from the short-distance part of 
the amplitude by expanding the short-distance part as a Taylor series in p. Second, one 
must decouple the Dirac indices and color indices that connect the short- distance part to 
the long-distance part. This can be accomplished by making use of Fierz rearrangements. 
In the factored decay rate, the long-distance parts correspond to the matrix elements of the 
NRQCD 4-fermion operators in the quarkonium state; the short-distance parts correspond 
to the imaginary parts of the coefficients of those operators in the NRQCD lagrangian. 

In order to see in more detail how the perturbative analysis leads to the results that we 
have obtained from NRQCD, let us consider two examples: annihilation of S-wave and P- 
wave quarkonium at leading nontrivial order in v and through order a^ in QCD perturbation 
theory. We use the specific example of decays into 2 and 3 gluons in the discussions below. 
However, the essential ingredients of the discussion apply also to decays into a light quark- 
antiquark pair and decays into a qq pair and a gluon. 



57 



B. Annihilation of S-wave Quarkonium 

The first step in the analysis of annihilation of S-wave quarkonium is to identify the short- 
distance part in the topological factorization of the amplitude. The dominant component 
of the bound-state wavefunction consists of a heavy quark and antiquark in a color-singlet 
state. We take the QQ pair to have total 4-momentum P and relative 4-momentum 2p. We 
assume that, owing to the bound-state dynamics, the Q and Q have inverse propagators 
(|P ± p)^ — M^ of order M^f ^, with v"^ <^ 1. In the meson rest frame, the energies \Pq ±po 
of the Q and Q then differ from the mass M by order Mv"^ and their momenta |P ± p are of 
order Mv. At order a^, the QQ pair can annihilate into two gluons through the diagrams in 
Fig. ^. By energy conservation, the two gluons must both have momenta of order M. At this 
order, the topological factorization of the annihilation rate is trivial. The QQ annihilation 
amplitude belongs entirely to the short-distance part of the annihilation rate in Fig. H, while 
the quarkonium wavefunction belongs to the long-distance part. 

We next consider the annihilation rate of the QQ pair at order a^. This rate has con- 
tributions from the annihilation into three gluons through the diagrams in Fig. ^ and also 
from the annihilation into two gluons, due to the interference between next-to-leading or- 
der diagrams such as those in Fig. ^ and the leading-order diagrams in Fig. ^. We begin 
by examining the infrared divergences in the diagrams for the emission of a real gluon of 
4-momentum / shown in Fig. ^. As we have already explained, the diagram in Fig. §(c) 
contains no infrared divergence. For the diagrams in Figs. §(a) and (b), we identify the 
infrared contribution that is leading in v by assuming that Pq ~ 2M, that /q = |1|) P, and 
p are of order Mv, and that po is of order Mv"^ . The emission vertex for the gluon with 
momentum I and the two adjacent heavy-quark propagators can then be approximated as 
follows: 



. ^A« 



{±\P + pf -M^ + ie ' (±iP + pTiy-M^ + 
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M(l±7o) „ M(l±7o) 
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[±\P + pY - M^ + ie ' -2Mlo + ie 
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M(l±7o) f ±< \ .... 

where the upper and lower signs correspond to Fig. ^(a) and Fig. ||(b), respectively. (In 
the case of the lower sign, the order of the gamma matrices should actually be reversed, 
but the last line is unaffected.) The factor ±g'^^ is called the "eikonal vertex", and the 
factor l/(— /o + ^e) = l/(— 11|) is called the "eikonal propagator". Their product is called the 
"eikonal factor" . We see that the eikonal factor for the contribution of Fig. §(a) is equal and 
opposite in sign to the eikonal factor for the contribution of Fig. @(b). All other propagator 
and vertex factors in the two diagrams are the same. If the QQ pair is in a color-singlet state, 
then the color factors in the two diagrams are also the same, and the infrared contributions 
from the region |1| ^ cancel. This cancellation is a consequence of the fact that, in the 
infrared limit, the soft gluon couples to the color charges of the quark and antiquark. Since 
the quarkonium is a color singlet, the quark and antiquark have opposite color charges. 

Because of the infrared cancellation, the topological factorization of the real-emission 
diagrams in Fig. ^ is trivial. The amplitudes for QQ -^ ggg all belong to the short distance 
part of the annihilation rate in Fig. 0, while the quarkonium wavefunction belongs to the 
long-distance part. 

Now let us turn to the virtual-gluon-emission diagrams shown in Fig. ^ Once again, we 
can identify the infrared part by neglecting / and p compared to M. As in the preceding 
example, the eikonal vertices are proportional to (?^° times the quark (or antiquark) charge. 
For Fig. |^(a), the eikonal propagator factors associated with the exchange of the gluon with 
momentum / are [l/(— /q + *e)] [—1/(^0 + ^e)]- Each of the diagrams of Figs, ^(b) and (c) 
contains a mass renormalization, which we subtract. The remaining contribution is a wave- 
function renormalization, half of which we absorb into the quarkonium wavefunction. The 
other half yields the eikonal propagator factor (— l/2)[l/(— /q + ie)]^- (The squared prop- 
agator appears after subtraction of the mass- renormalization contribution.) The eikonal 
factors from the three diagrams would cancel, were it not for the ie's in the propagator 
denominators. Instead, the eikonal factors yield 
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—, -1 ■■ ] -—, = 2m5{k)— -. (5.2) 

—to + It Iq + It —Lq + It —to + It — io + ^e 

The (5(/o) contribution arises because of the pinch in the Iq integration contour in the con- 
tribution of the diagram of Fig. ^(a). This 5{lo) contribution corresponds to the exchange 
of a space-hke gluon with temporal polarization between the quark and antiquark. That is, 
it corresponds to the Coulomb scattering of the quark and antiquark. Note that the factor 
multiplying 5{Iq) is divergent at Iq = 0. This somewhat unexpected divergence has arisen 
because we have neglected the relative momentum p of the heavy quark and antiquark. 
Had we retained that momentum in the propagator denominators, we would have obtained 
a 1/f singularity, where v is the relative velocity of the quark and antiquark. [This 1/v 
contribution is calculated in detail in (|A20|) -( |A22|) .] Ordinarily, in the absence of a collinear 



singularity, the phase space for two partons to be moving with small relative velocity would 
be unimportant. Here, that region of phase space is important by virtue of the quarkonium 
bound state. (In fact, it is the 1/v singularity that builds up the bound-state wave function 
in a perturbative analysis of the Bethe-Salpeter equation for positronium.) 

At this point, the topological factorization of the virtual-emission diagrams can be carried 
out. For the diagrams in Fig. |, one factors the following contributions into the long- 
distance part: the wavefunctions, with which we associate the square root of each quark 
or antiquark wavefunction renormalization, and the 1/v singularity that arises from the 
diagram of Fig. |^(a). The remaining contributions from these diagrams are factored into 
the short-distance part. 

Many discussions of perturbative factorization make use of the Grammer-Yennie tech- 
nique [|22| for analyzing infrared divergences. As an aside, let us indicate briefly how that 



technique would apply to the example at hand. From our previous discussion, we see that, 
in the infrared limit, the infrared-gluon vertex V^ is well approximated by q'^^Vq. Then we 
can write 

yM ~ gf^o XlL^ (5.3) 

k-it 

provided that /q is not small compared with the other components of /. This is always 

60 



the case for real emission. For virtual emission, we can eliminate the region of small /q by 
deforming the Iq contour of integration into the lower half of the complex plane. As can 
be seen from an examination of the propagator denominators, all of the singularities in the 
lower half of the complex plane are order M° away from the origin, except in the diagram 
of Fig. ^(a). This is the momentum-space manifestation of the fact that a space-like gluon 
(with small Iq) can be exchanged causally only between co-moving particles. In carrying out 
the contour deformation for the diagram of Fig. |^(a), and only in this case, we unavoidably 
pick up the residue of a pole at /q ~ 0. This residue yields the 1/v singularity in (|5.2| ) that 
was noted earlier. Along the deformed contours, the Grammer-Yennie approximation ( ^.3|) 
is valid. Substituting ( ^.3|) for the infrared vertices, we can make use of Ward identities 
(current conservation) to show that the contributions of the deformed contours cancel. In 
order to obtain this Ward-identity cancellation, one needs, in addition to the Grammer- 
Yennie contributions of Figs. ^ and ^ Grammer-Yennie contributions in which the infrared 
gluon attaches to the short- distance part of the process. But, as we have already argued, 
these diagrams give contributions that vanish in the infrared region, so there is no harm in 
applying the Grammer-Yennie approximation to them. 

After topological factorization, the short-distance and long-distance parts of the anni- 
hilation rate are still tied together by integrations over the relative 4-momenta p and p' of 
the QQ pairs entering and leaving the short distance part and by sums over the color and 
Dirac indices associated with the heavy-quark propagators. To complete the factorization, 
we must decouple these integrals and sums. 

The decoupling of the integration over p and p' is accomplished simply by expanding the 
short- distance contribution in a Taylor series in p and p'. Taking p and p' to be of order Mv 
and po and Pq to be of order Mv'^, we see that the Taylor expansion of the short- distance 
part corresponds to an expansion in powers of v. All of the dependence on p and p' is now 
in the long-distance part and in the explicit powers of p and p' from the Taylor expansion. 
To analyze S-wave decays at leading order in f , we need keep only the zeroth order terms 
in the Taylor expansion, which amounts to setting p = p' = 0. In the meson rest frame, the 
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QQ pair has total energy Pq which differs from 2M by an amount of order Mf ^ , and total 
momentum P of order Mv. At leading order in f , we can set Pq = 2M and P = in the 
short- distance part of the annihilation rate. Thus the incoming quark and antiquark can be 
taken to be on their mass shells and at threshold. 

The decoupling of the color indices connecting the short-distance and long-distance parts 
of the annihilation rate is straightforward. The short-distance part is a color tensor Cij^ki, 
with color indices i and j for the incoming Q and Q and j and k for the outgoing Q and Q. 
The indices i and j can be decoupled from the tensor by using the Fierz rearrangement 



1 



^i'i ^jj' = W ^Ji ^i'J' + 2 T^i T^j" (5.4) 



A similar rearrangement can be used for the indices k and /. By color symmetry, T°':Cij,kk 
must vanish and T°',Cij^kiTij. must be proportional to the unit tensor 5°'^. The resulting 
rearrangement formula is 

CijM = j^ Sji {Ci'i'j'j') 5ki + j^2 _i ^ji [Ti'j'Ci'j',kn'Ti,k,j T^i. (5.5) 

The indices have been decoupled from the tensor by decomposing it into a term in which 
both pairs of indices are projected onto a color-singlet state and a term in which both pairs 
of indices are projected onto a color-octet state. For S-wave quarkonium, the dominant 
Fock state contains a color-singlet QQ pair, so only the first term on the right side of ( ^.5|) 
contributes at leading order in v. 

The Dirac indices connecting the short- distance and long-distance parts of the amplitude 
can be decoupled in a similar way, although the algebra is a little more cumbersome than 
it is for the color indices. Having set p = p' = and P = (2M, 0) in the short- distance part 
of the amplitude, we find that the numerators of the four propagators connecting it to the 
long-distance part reduce to M(7o + 1) for the quarks and M(— 70 + 1) for the antiquarks. 
The Dirac structure of the short- distance part of the amplitude is therefore a tensor Tij^kh iii 
which the Dirac indices i and k of the quarks are contracted with projectors P+ = (l + 7o)/2, 
and the Dirac indices j and / of the antiquarks are contracted with projectors P_ = (1— 7o)/2. 
The indices i and j can be decoupled from the Dirac tensor by using the Fierz rearrangement 
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(P+U iP-)n' = \ (75^),, (75P-),,v + \ (a^5P+),, (^"75^-),,,. , (5.6) 

where a* = (i/2)e°*-''^[7'', 7^^]. A similar Fierz rearrangement can be used to decouple the 
indices k and / from the Dirac tensor. Since all 3-momenta have been set to in the Dirac 
tensor F*-^''^', there is no 3-vector on which F*-'''^' can depend. Its transformation properties 
under rotations then imply that the vector {a°''^^P-)ijVij^ki{lbP+)ik must vanish, while the 
tensor (o"'^75-P_)ijFjj^fci(a^75P+)ifc must be proportional to the Cartesian unit tensor 5"*. 
Consequently, one obtains the rearrangement formula 

{P+)i'i iP-)jj' ^i'j',k'l' {P+)kk' iP-)l'l 

= -^ il5P+)ji [i-f5P-)i,jSi'j',k'l'il5P+)i,k] (75^-)h 

+ ^ (a"75P+),. [(^'^75^-),,, F,,, ,,, (^'^75^),,,] (^'"75P-)fc, • (5.7) 

This rearrangement of the Dirac indices corresponds to the decomposition of the Dirac tensor 
into spin-singlet and spin-triplet pieces. The Dirac matrix 75P- = -P+75 projects a QQ pair 
at rest onto a state with total spin 0, as can be seen from the identity 

^(0,0|i,m;i,m') M„t;^/ = — = (75P_) , (5.8) 

mm' V ^ 

where Um and Vm' are Dirac spinors evaluated at zero 3-momentum. Similarly, the Dirac 
matrix cr*75P_ = P^a^'-y^ projects a QQ pair at rest onto a state of total spin 1: 

Y,{l,M\\,m-\,m')u^Vm' = -j=U^'' {cr'j,P_) , (5.9) 

mm' V -^ 

where f/*^* is the unitary matrix that transforms from the Cartesian basis to the basis of 
angular-momentum eigenstates. 

Now that we have decoupled the integrations over p and p' and the sums over color and 
Dirac indices, the factorization of the annihilation rate is complete. In the rearrangement 
identity ( |5.7| ), the factors on the right side that are enclosed in square brackets belong to the 
short-distance part of the annihilation rate. They correspond to the operator coefficients 
in the NRQCD approach. The remaining factors to the right and to the left of the square 
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brackets belong to the long-distance part. It is evident that the long-distance parts of the 
annihilation rate can be reproduced by matrix elements of local operators in the quarkonium 
state. The two operators that contribute to the annihilation of S-wave states at leading order 
in V may be identified as 

-^75P+^^75P-^ ~ Oi{^Sq), (5.10a) 

-^o-75P+^-^o-75P_^ ^ Ci(35i), (5.10b) 

where \1/ is the Dirac field for the heavy quark. For matrix elements between quarko- 
nium states, these operators reduce at leading order in v to the NRQCD operators 
C^i(^'S'o) = V^^XxV and 01(^51) = Tp'^crx ■ X^^^i^^ respectively. Thus, perturbative factoriza- 
tion yields the same operator matrix elements as appear in the NRQCD analysis. It should 
be noted, however, that the identifications (|5.10|) are not unique. For example, the operator 
— ^75\I'\1'75\E', when sandwiched between quarkonium states, also reduces at leading order 
in V to 01(^5*0) and both -'^ct'^^^ ■ ^a'^^^ and -"^7^ ■ '^7\I^ reduce to 01(^5*1). 

C. Annihilation of P-wave Quarkonium 

Now let us analyze the annihilation of P-wave quarkonium at leading nontrivial order 
in V. First we note that, because the spatial part of the P-wave quarkonium wavefunction 
transforms under rotations like a vector, the p-independent part of the QQ annihilation 
amplitude vanishes on carrying out the angular part of the integration over p. Thus, we 
must retain terms with one factor of p in the annihilation amplitude, which means that the 
leading amplitude is down by one power of v relative to the S-wave case. 

At order a^, the annihilation proceeds through the diagrams in Fig. ^. In this case, the 
factor of p in the QQ annihilation amplitude can come only from expanding the propagator 
of the virtual heavy quark, which is off its mass-shell by an amount of order M. The 
topological factorization is therefore trivial. The amplitude for QQ -^ gg belongs to the 
short-distance part of Fig. 0, and the quarkonium wavefunction belongs to the long-distance 
part. 
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We next consider the annihilation at order af, which receives contributions from the 
real-emission diagrams in Fig. ^ and from the virtual-emission diagrams in Fig. ^. The 
factor of p can come from one of two sources: the purely short-distance (infrared-safe) part 
of the diagram, or the potentially infrared- divergent part, which consists of the soft gluon 
and the heavy-quark propagators to which it attaches. 

If the factor of p comes from the short- distance part of the diagram, then the analysis 
of the infrared divergences goes through exactly as in the S-wave case. Infrared divergences 
cancel between the real-emission diagrams, but the exchange of a virtual gluon between the 
Q and Q [Fig. 0(c)] results in a 1/v singularity. Topological factorization is trivial, except 
for this 1/v singularity. It must be factored into the long-distance part of the annihilation 
rate. 

We proceed to consider the case in which the factor of p comes from the potentially 
infrared-divergent part of the diagram. We consider separately the cases of virtual-gluon 
emission and real-gluon emission. 

The diagrams for virtual-gluon emission are shown in Fig. 0. The potentially infrared- 
divergent part of the amplitude includes the factors in the first line of ( [5.1| ) . The required 
factor of p can come either from a p ■ 7 in the numerator of a propagator or from ex- 
panding out the denominator. The terms with a factor of p that comes from a propagator 
denominator are easily seen to be suppressed by a power of v. The terms that contain a 
p ■ 7 in the numerator are also suppressed by a factor of v because of the Dirac structure. 
To see this, first consider the case of a soft gluon with temporal polarization. From the 
identity P+p ■ 7 = P ■ 7-P-, one can see that the factor p • 7 connects "large" components 
of Dirac matrices to "small" components. This gives rise to the suppression by a factor of 
V. Now consider the case of a virtual gluon with spatial polarization vector e. Both of the 
spatial-gluon vertices bring in factors of e ■ 7. The combined effect of these two factors and 
the factor of p • 7 is again to connect large and small components, which costs a factor of 
V. Thus for virtual-gluon emission, there are no infrared divergences at leading order in 
V. The topological factorization is therefore trivial. The amplitude for QQ — > gg belongs 
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entirely to the short- distance factor in the annihilation rate, aside from the square root of 
each heavy-quark or heavy- antiquark wavefunction renormalization, which we associate with 
the quarkonium wavefunction. 

Finally, we consider the case of real-gluon emission through the diagrams in Fig. ^ Let 
us examine the infrared limit of the diagrams in Fig. ^(a) and ^(b) for the case in which the 
soft gluon with momentum / has spatial polarization vector e. The emission vertex and the 
adjacent heavy-quark propagators can be approximated as follows: 
(±ip + p)-7 + M (±|P + pT/)-7 + M 



e-7 



±iP + p)2-M2 + ie ' {±^P + pTiy-M^ + ie 



Af (1 ± 7o) - p ■ 7 M(l ± 7o) - P ■ 7 

e -7 



(±|P + p)2 - M2 + te ' -2M/o + ie 
M(l±7o) /2p-e\ 



(5.11) 



(±|P + p)2-M2 + ie \^ -|1| ; ■ 
The upper and lower signs apply to Figs, ^(a) and^(b), respectively. (For the lower sign, the 
order of the Dirac matrices should actually be reversed.) In the last line of (|5.11 ), we have 



retained only those numerator terms that contain one power of p. The factor 2p • e/( — 11|) 
is the infrared-emission factor. In contrast with the S-wave case, the infrared contributions 
from the two real-emission diagrams add, rather than cancelling. Because we have retained 
one power of p, the soft gluon couples to the color current of the heavy quark, rather than 
to the color charge. Since the heavy quark and antiquark have opposite color charges and, 
in the CM frame, opposite momenta, their color currents are equal. Note that, because 
of the vector p in the infrared-emission factor, the emission of the soft gluon changes the 
orbital-angular- momentum quantum number of the QQ pair by one unit, but it does not 
flip the spin of the quark or antiquark. Thus, it converts the heavy quark and antiquark 
from a color-singlet P-wave state to a color-octet S-wave state. 

In the decay rate, we must integrate the infrared emission factors from the square of the 
sum of the amplitudes over the phase space of the gluon. Keeping only the logarithmically 
divergent part of the integral, we flnd the result 
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where A is an infrared cutoff of order Mv, and we have arbitrarily set the upper hniit on |1| 
to M. We have introduced a factorization scale A to separate the infrared divergence from 
the short-distance part of the integral. The long-distance contribution that is proportional 



to log(A/A) in ( |5.12| ) can be interpreted as the probability for a heavy quark and antiquark 
in a color-singlet P-wave state to make a transition to a color-octet S-wave state by radiating 
a soft gluon. 

We can now carry out the topological factorization of the diagrams in Fig. § for real gluon 
emission. The square of the amplitude for QQ — > ggg, integrated over phase space, belongs 
to the short-distance part of the annihilation rate in Fig. |^, except for the second term 
on the right side of ( |5.12|) . This term, which contains the infrared divergent contribution 



that arises from the emission of the soft gluon in Figs, ^a) and §(b), is included in the 
long-distance part, along with the quarkonium wavefunction. The soft gluon is an example 
of a light parton that connects the initial and final wavefunctions in Fig. |[ Note that, in 
this contribution to the annihilation rate, the heavy quark and antiquark enter the short- 
distance part in a color-octet S-wave state. We call this contribution to the annihilation rate 
the "color-octet contribution" . In all the other contributions to the P-wave annihilation rate 
at this order, the heavy quark and antiquark enter the short- distance part in a color-singlet 
P-wave state. We refer to those contributions as the "color-singlet contribution" . 

At this point, we have identified the long- and short-distance parts in the topological 
factorization of the annihilation rate. It remains only to decouple the integrations over the 
relative momenta p and p' of the QQ pairs and the sums over color and Dirac indices. 

We first discuss the color-singlet contribution. The color indices of the short- distance 
and long-distance parts are easily decoupled by using the rearrangement identity (^.51). Only 
the first term on the right side of ( |5.5| ) contributes, since the QQ pair is in a color-singlet 
state. In order to decouple the integrals over the relative momenta p and p' of the initial 
and final QQ pairs, we expand the short distance part as a Taylor series in p and p' . At 
leading order in f , we set po = Po' = and keep only those terms linear in both p and p'. 
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The resulting amplitude has the structure p^'^YjkiP' ■> i^ which the Dirac indices i and k of 
the quark are contracted with projectors P+ and the Dirac indices j and / of the antiquarks 
are contracted with projectors P_. The Dirac indices of the initial and final QQ pair are 
decoupled from the short distance factor F^"^; by applying the Fierz identity ( [5 .61 ) to both the 
initial and final indices. The resulting rearrangement identity can be greatly simplified by 
making use of the fact that the Fierz-decoupled short-distance parts transform like tensors 
under rotations and the fact that there are no three-vectors on which V^^i can depend. For 
example, the tensor (75P+)jjT^'^;(75P+)/fc must be proportional to the unit tensor 5"*", and 
the tensor {cr°''~^^P^) ijV^'^f,i{^^P^)ii^ must be proportional to the totally antisymmetric tensor 
^amn ^ Siuce (cr"75P+) jjr^^'^;(o"^75P+);fc is a Cartesian tensor in 3 dimensions with 4 indices, 
it can be decomposed into a linear combination of the three tensors (5"™'^''", ("-"^^^^^^ ^ and 
ij^^afe^mn _|_ ^an^mfe-j _ i^am^fen^ which corrcspoud to total augular momentum 0, 1, and 2, 
respectively. Consequently, one obtains the rearrangement formula 
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-- ((p X ^)-75P+),. [e'^"^ K75P-).,, F^f^., ,,, (75P+),,,] (P'™75P-),, 



24 
+ 1 (p™75P+),. [e'^"^ (75P-),,v F^,^,,,,, (a^75n)..'] ((P x (t)'"75P-)., • (5.13) 

We use the notation T^""^^ for the symmetric traceless part of a tensor T"* and T''*''] = 
|(T'*'' — T^"') for the antisymmetric part. 

At this point, the factorization of the color-singlet contribution to the annihilation rate 
is complete. The factors in square brackets on the right side of ( |5.13| ) belong to the short- 
distance part, while the factors to the right and to the left of the square brackets belong 
to the long-distance part. It is apparent from the rearrangement identity ( ^.13| ) that the 
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long-distance parts can be reproduced by matrix elements of the following local operators: 

-^(-fD)75P+^-^(-fD)75P-^ ^ Oi('Pi), (5.14a) 

-|^(-fD-o-)75P+^^(-fD-o-)75P_^ ~ Oi('Po), (5.14b) 

-i^(-f D X ^)75P+^ ■ ^(-fD X o-)75P_^ ^ Oi('Pi), (5.14c) 

-^(-fl3(V^))75P+^^(-fB(V^))75P-^ ~ O1CP2), (5.14d) 

-^(-f D)75P+^ • ^(-f D X o-)75P-^, (5.14e) 

-^(-f D X o-)75P-^ ■ ^(-f D)75P+^. (5.14f) 

The matrix elements of the last two operators vanish for a quarkonium state that is a charge- 
conjugation eigenstate. The other four operators reduce at leading order in v to the NRQCD 
operators C(^Pi) and 0{^Pj), J = 0, 1, 2, respectively. 

Finally, we consider the factorization of the color-octet contribution to the annihilation 
rate, for which the short- distance part involves the annihilation of a QQ pair in a color- 
octet S-wave state. The color indices of the short-distance and long-distance parts are easily 
decoupled by using the rearrangement identity ( p.5[ ). Only the second term on the right side 
of ( p. 51 ) is non-vanishing for the color-octet contribution. The decoupling of the momentum 
integrations and the Dirac indices proceeds along the same lines as for S-wave quarkonium, 
which was discussed in subsection [VB| . The momentum integrations are decoupled by 
Taylor-expanding the short- distance part in p and p', and setting p = p' = 0. The decoupling 
of the Dirac indices is accomplished by using the rearrangement formula ( ^.7|) . The factors in 
square brackets in (|5.7| ) belong to the short- distance part of the annihilation rate, while the 
factors to the right and to the left belong to the long-distance part. From the rearrangement 



identities ( |5.5| ) and (|5.7|) , it is evident that the long-distance parts are reproduced by matrix 
elements of the operators 

-^75T"P+^^75T"P_^ ^ OsCSo), (5.15a) 

-^<T75T"P+^ ■ ^cr75T"P_^ ^ OsCSi). (5.15b) 

These operators reduce, at leading order in v, to the operators Os(}So) = ip'^T'^xX^T'^i^ ^"^^ 
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C^si^Si) = iP^(tT"'x ■ X^(^T°"ip in the NRQCD analysis. The matrix elements of Os(^So) and 
Osi^Si) include the probability factor proportional to log(A/A) in (|5.12|). The logarithmic 



dependence on A is reflected in the evolution of these operators, which is given in ( p. 28a] ) 
and ( |3.28b| ). Thus, the factorization scale A in the perturbative approach can be identifled 
with the ultraviolet cutoff of NRQCD. 
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VI. PRODUCTION OF HEAVY QUARKONIUM 

In this section, we present a general factorization formula for computing inclusive heavy- 
quarkonium production rates in high-energy processes that involve a momentum transfer 
Q^ that is of order M^ or larger. In the case of S-wave quarkonium, our factorization 
formalism coincides with the "color- singlet model" for quarkonium production |^ in the 
nonrelativistic limit, but it also allows the systematic calculation of relativistic corrections 
that are suppressed by powers oi v. In the case of P-wave quarkonium, our formalism 
reveals that the color-singlet model is incomplete, even at leading order in f , and must be 
supplemented by including the "color-octet mechanism" for P-wave quarkonium production 



A. Factorization of the Production Rate 

Our goal, as in the discussion of heavy-quarkonium annihilation, is to express the inclu- 
sive production rate for a quarkonium state in a factored form. That is, we wish to write 
the production rate as a sum of terms, each of which consists of a short-distance part, which 
can be calculated in QCD perturbation theory, multiplied by a long-distance part that can 
be expressed as a matrix element in NRQCD. Our arguments for the factorization of the 
production rate are based on the all-orders properties of QCD perturbation theory. In this 
sense, the level of rigor of these arguments is comparable to that in the proofs of factoriza- 



tion for the Drell-Yan process for lepton pair-production in hadron-hadron collisions p5 
These arguments are less rigorous than those that we have given for the factorization of 
the quarkonium annihilation rate. The latter arguments rely only on the general space- 
time structure of the annihilation process and on the validity of the effective-field-theory 
approach. Their level of rigor is comparable to that in the proofs of factorization in deep- 
inelastic lepton-hadron scattering, which can be formulated in terms of the operator-product 
expansion. 
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When a quarkonium state is produced in a process that involves momentum transfer 
Q^ of order M^ or larger, the production of the QQ pair that forms the bound state takes 
place at short distances of order \jM or smaller. A simple example of such a process, 
which the reader can keep in mind throughout the following discussion, is the production in 
e^e~ annihilation at a center-of-mass energy y^ ^ M of a heavy quarkonium iJ, with 4- 
momentum P, recoiling against two light hadron jets. At leading order in QCD perturbation 
theory, the relevant parton process is e^e" -^ QQgg- We take the Q and Q to have momenta 
P/2 + p and P/2 — p. The relative 3- momentum p must be of order Mv in the P = frame 
in order for the QQ pair to have a significant probability for forming the bound state H. 
The amplitude for the production of the QQ pair is insensitive to changes in the relative 
4-momentum p that are much less than M, and therefore the quark and antiquark are 
produced with a separation of order 1/M or less. Similarly, the square of the amplitude 
is insensitive to changes in the the total 4-momentum P of the heavy pair that are much 
less than M. Thus, the product of one amplitude and the complex conjugate of a second 
will contribute significantly to the QQ-production cross section only if the corresponding 
production points are separated by a distance of order 1/M or less. We therefore conclude 
that the production of the QQ pair is indeed a short-distance process that takes place within 
a distance of order 1/M. 

In the framework of NRQCD, the effect of the short- distance part of a production ampli- 
tude is simply to create a QQ pair at a spacetime point. The formation of the quarkonium 
state H from the QQ pair takes place over distances that are of order l/{Mv) or larger in the 
quarkonium rest frame, so it is described accurately by NRQCD. Therefore, in NRQCD, the 
production rate (the square of the amplitude summed over final states) involves the creation 
of a QQ pair at a spacetime point, its propagation into the asymptotic future, where the 
out state includes the quarkonium H, and, finally, the propagation of the QQ pair back in 
time to the creation point. That is, the long-distance part of the production rate is given 
in NRQCD by vacuum matrix elements of local 4-fermion operators. The effects of the 
short-distance parts of the production rate are taken into account through the coefficients 
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of the 4-ferinion operators. Since the final state must include a quarkonium, the 4-ferniion 
operators that appear in production cross sections involve projections onto the space of 
states that contain, in the asymptotic future, the quarkonium state H plus anything else. 
The generic form of a production operator is 

Of = x^/c^v^ f E E 1^ + X) {H + X\\ v^t/cix 

\ X mj / 

= X^K,n^(a^HO'H)^'^K:'^X, (6.1) 

where the sums are over the 2J + 1 spin states of the quarkonium H and over all other final- 
state particles X. In the second line of ( |6.1| ), the projection has been expressed compactly 
in terms of the operator a\j that creates the quarkonium H in the out state. A sum over 
the angular- momentum quantum numbers mj is implicit in a/^an- The factors /C„ and /C^ 
in the operator are products of a color matrix (either the unit matrix or T*^), a spin matrix 
(either the unit matrix or cr*), and a polynomial in the covariant derivative D and other 
fields. The overall operator O^ is invariant under color and spatial rotations.a We assume 
that any matrix elements of O^ will be evaluated in the quarkonium rest frame; otherwise 
the factors /C^ and /C^ may depend on the 4-momentum of the quarkonium. 

It is convenient to introduce notation for the production operators that is analogous 
to that for the decay operators defined in ( |2.10| ) and ( p.l2[ ). The production operators of 
dimension 6 are 

C?f('5o) = xV(aW)V'^X, (6.2a) 

Of (^Si) = x^V (a^^an) V^V^x, (6.2b) 

O^i^So) = X^T> (a^H^H) ^^T'^X, (6.2c) 



^Here we consider explicitly only unpolarized production of heavy quarkonium. In the case of 
polarized production, a\j would create a state of definite polarization, and Kn and /C^ would, 
in general, depend on one or more vectors associated with the incoming particles, such as the 
directions of their spins and momenta. 
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Oi{'S^) = x^cr'T'^ij {a^^an) ^^a^T'^x- (6.2d) 
Some of the color-singlet production operators of dimension 8 are 

OU'Pi) = X\-01^ {aW) i^H-0')x, (6.3a) 

O^'Po) = I X^-fD ■ a)^ (aW) ^^(-fD ■ ct)x, (6.3b) 

O^CPi) = \ X^(-tD X a)V (aW) ^^(-f D X aYx, (6.3c) 

OU'P2) = xT(-f B(V^))^ (a^H<^H) ^^(-1 B(V^))x, (6.3d) 

V^i'S,) = ^ [xV (^a^) i^K-^^fx + h-c] , (6.3e) 

^f e^i) = ^ [x^^V (4«h) ^V^(-|D)2x + h.c] . (6.3f) 



Given that the long-distance part of the production rate can be expressed in terms of 
vacuum matrix elements of operators of the form given in (|6.1| ), the inclusive production 
cross section must have the form 

-(H) =E§^(0|a^(A)|0), (6.4) 

n 

where it is understood that the matrix element is to be evaluated in the quarkonium rest 
frame. The short-distance coefficients F„ depend on all the kinematic variables of the pro- 
duction process, but they are independent of the quarkonium state H. Equation ( |6.4|) is the 
equivalent for production of our factorization formula ( p.l4| ) for quarkonium decay. 



Beyond leading order in perturbation theory, interactions involving soft (infrared) gluons 
and gluons coUinear to the final-state jets potentially spoil this factorization picture, both 
by making the QQ-production process long-ranged and by making connections between the 
outgoing quarkonium and the final-state jets that destroy the topological factorization. In 
the case of quarkonium decay, we were able to use the KLN theorem to argue that such 
final-state soft and collinear interactions cancel in the inclusive decay rate. In the case of 
quarkonium production, the KLN theorem does not apply directly because we have specified 
that the final state contain the quarkonium: some of the cuts in the KLN sum are missing. 
Cuts are missing only for diagrams in which a soft or collinear gluon attaches to one of the 
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heavy Q or Q lines. If only one end of a gluon attaches to a Q or Q line and the other end 
attaches to a final-state jet, then the sum over cuts along the jet line is sufficient by itself 
to effect the KLN cancellation. If both ends of a soft gluon attach to a heavy Q or Q line, 
then there is no KLN cancellation. However, this contribution is part of the matrix element 
of the NRQCD 4-fermion operator. 

In the case that O^ is a color-octet operator, one might worry that, because the inter- 
mediate states in the first line of (|6.1|) carry net color charge, the factorization of the cross 



section in ( |6.4|) is not valid. Owing to the property of confinement, such colored states have 
infinite energy. (Their energies would be finite in a finite volume, however.) Of course, 
the complete final state is color neutral and contains only color-singlet hadrons. One can 
picture the color neutralization of the partons in perturbation theory as a process involving 
soft-gluon exchanges between the partons. In particular, there can be color-neutralizing 
soft-gluon exchanges between partons that are comoving with the quarkonium and partons 
in other hadron jets produced by the short distance process. However, the KLN argument 
tells us that, at least in perturbation theory, the infrared and collinear divergences from such 
soft interactions cancel in the inclusive quarkonium production rate. That is, for purposes 
of computing the inclusive quarkonium production rate, the colored partons can be treated 
as if they were unconfined. Of course, the complete operator O^ is invariant under color 
rotations, and one can deal with it without referring to the troublesome colored intermediate 
states by making use of the form given in the second line of (|6.1| ). This approach might be 
useful in lattice measurements of the production matrix elements. 

If we consider production of quarkonium in hadron-induced processes, then a host of 
new difficulties arise in proving that the production rate factors. These include exchanges of 
soft, collinear, and Glauber (quasi-elastic) gluons involving spectator partons in the initial 
state and exchanges of soft and collinear gluons involving active partons in the initial state. 
Rather than discuss the resolution of these difficulties here, we will merely assume that 
the Glauber divergences cancel, that the only noncancelling infrared divergences are those 
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associated with the matrix elements of the 4-fermion operators, and that the noncancelhng 
colhnear divergences can be absorbed into initial parton distributions. We refer the reader 
to the proofs of factorization of the Drell-Yan cross section p5| , |26[| for detailed discussions 
of these points. Given these assumptions, the factored form ( |6.4| ) holds to all orders in 
perturbation theory. It should be noted that, in the case of hadron-hadron collisions, there 
is a limit to the precision of the factored form of the cross section. Generally, because of 
soft exchanges between spectators, one can prove only that a factored form holds through 
next-to-leading order in an expansion in inverse powers of the large momentum transfer Q^ 



|2q| . Beyond that order, factorization is known to fail ||27| 



B. Relation of Production Matrix Elements to Decay Matrix Elements 



The NRQCD matrix elements that appear in the production rate ( |6.4|) are related to the 
NRQCD matrix elements that appear in decay rates through a crossing of the quarkonium 
from the final state to the initial state. This relation is analogous to the one between parton 



distribution functions and parton fragmentation functions [^. In general, the crossing 
relation is very complicated. There are, however, two instances in which one can obtain 
simple results. 

Through order a^ in QCD perturbation theory, the crossing relation between (i7|0„|i7) 
and the corresponding production operator (0|(9^|0) is a simple equality, up to a factor of 
2 J + 1 for the number of spin states. Finite-order perturbation theory is usually of little help 
in dealing with long-distance matrix elements. It does tell us, though, that, to leading order 
in Qfs, the evolution equations for the production operators are the same as the evolution 
equations for the corresponding decay operators. For example, the evolution equation for the 
production matrix element (0|C8''(^5'o)|0) in terms of (0|(9i''(^Pi)|0) is identical at leading 
order in a^ and in v to the evolution equation (p.28a|) for the corresponding decay matrix 
elements: 

A^(0|08'^C5o)|0) = ^]^(0|Of^CP0|0). (6.5) 
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When O^ is a color-singlet operator, the vacuum-saturation approximation can some- 
times be used to simplify the matrix element. Assuming that the sum over states in the first 
line of (|6.1|) is dominated by the quarkonium state H plus the vacuum, we obtain 



^i2J+l){H\0„\H), (6.6) 

where C„ = ip"^ K.'nXx'' ^n'ip ■ In the second line, we have used the rotational invariance of the 
operator ip'' IC'^xl^) i^lx'^ ^n'ip , which implies that the matrix element is identical for each of 
the 2 J + 1 angular-momentum states H that differ only in the quantum number rrij. In 
the last line, we have used the vacuum-saturation approximation (p.8|) for the decay matrix 
element {H\On\H). 

For the vacuum-saturation approximation to be a controlled approximation, we must be 
able to show that the contributions of all the other states in the sum in ( |6.1| ) are suppressed 
by powers of v. This is in fact the case if the operator (9^ creates and annihilates the QQ 
pair in the angular-momentum state that corresponds to the dominant Fock state of the 
meson H. In this case, the vacuum-saturation approximation result ( |6.6| ) is correct up to 
an error of relative order f ^. 

In the case of a color-octet operator, the states \H + X) in the first line of (|6.1|) have 
nonzero color, and the vacuum-saturation approximation is not applicable. In perturbation 
theory, we can approximate the sum by retaining only the terms involving intermediate 
states \H + g) that contain a single gluon. Similarly, we can approximate the sum for 
the corresponding decay matrix element by retaining the terms that involve single-gluon 
intermediate states \g). The resulting matrix elements {0\x^ JCn^PlH + g) and {g\x^ }Cnip\H) 
are related by crossing. Unfortunately the crossing relation is a simple equality only at 
leading order in perturbation theory. In the absence of any rigorous relation between them, 
we treat the matrix elements of the color-octet production operators and the color-octet 
decay operators as independent nonperturbative quantities. 
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C. Computation of the Operator Coefficients 

The short-distance part of the quarkonium production rate is insensitive to the long- 
distance QQ dynamics. Therefore, following the same reasoning as in Section |11H|, we can 



exploit the equivalence of perturbative QCD and perturbative NRQCD at long distances 
as a device to calculate the coefficients of the matrix elements in ( |6.4| ). We compute the 
production rate for an on-shell QQ pair with small relative momentum using perturbation 
theory in full QCD. Then we use perturbation theory in NRQCD to compute the matrix 
elements of 4-fermion operators O^^ , which are analogous to those in ( [Ol) except that the 
projection is onto on-shell QQ states. The short- distance coefficients are then determined 
by the matching condition 



a 



{QQ) 



pert. QCD n 



E#^(o|(^„«^(A)|o) 



(6.7) 

pert. NRQCD 



By expanding the left and right sides of (|6.71 ) as Taylor series in the relative momentum 
p between the Q and Q, we can identify the coefficients of the individual operators. They 
correspond to the infrared- and collinear-finite parts of cross sections for QQ production. 



One useful way to evaluate the left side of (|6.7] ) is to express the projection of the product 



u{P/2 + p)v{P/2 — p) of the Q and Q spinors onto a particular angular momentum state 
in Lorentz-invariant form. We refer the reader to Ref. ||2^ for examples. Then the left side 
of ( |6.7| ) can be evaluated in any convenient frame, such as the CM frame of the overall 
production process. It is understood, of course, that the matrix elements on the right side 
of ( |6.7|) are to be evaluated in the rest frame of the quarkonium. 

D. S-wave Production 

We now apply the factorization formalism to the production of S-wave quarkonium 
through relative order f^. For definiteness, we use the lowest-lying S-wave levels of charmo- 
nium for the purpose of illustration. Of course, the results that we give generalize imme- 
diately to other S-wave quarkonium systems. According to ( |6.4|) , the cross section for the 
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inclusive production of S-wave charmoniuni is 

^(^c) = ^^^ mOTCSoM + ^^ mvrCSom + 0{v'a\ (e.Sa) 

a(^) = ^^^ {^\Ot('S,m + ^^(0|Pf r5i)|0) + 0(t;V). (6.8b) 



The vacuum-saturation approximation ( |6.6|) can be used to reduce the 4-fermion matrix 
elements to products of matrix elements between the vacuum and the quarkonium state. 



These can, in turn, be related to the quarkonium wavefunctions given in Section [IIID . 



Finally, heavy-quark spin symmetry can be used to reduce the matrix elements to the same 
three nonperturbative parameters that appear in charmonium decay: |-R,yJ^, l-R^P; ^.nd 



Re(i?5'* V^-Rs). Taking into account factors of 2J + 1 for the number of spin states, we find 
that the cross sections are 

^(^c) = o^/i^2 ^vc n^.jA Re(i?5 Vms) + 0{v a), (6.9a) 



27rM2 



27rM4 



3iV,Fi(3^i) —2 3N,G^i^S,) ^ ,_ 



27rM2 



^(^) = L ^2 ^'A L:;4 MRs^'^'Rs) + 0{v'a). (6.9b) 



27rM4 



If we require only accuracy to leading order in v, then we can simplify the production 



rates in (|6.9|) further by dropping the terms proportional to Re(i?s* V^-R^) and replacing R^^ 
and R^ by their weighted average Rs- We then recover the familiar factorization formulas 
used in most previous work: 



NcFA'S, 



0) 



2 



^(ric) = ^'y^^^r' Rs ' + 0(t;V), (6.10a) 



3iV,Fi(^5i 



2 



<y{^) = "^'Hm^^^' ^' ' + ^^'''''^- ^^-^^^^ 



In applying the factorization formula ( |6.4D , one should keep in mind that the short- 
distance coefficients -F„(A) depend not only on as{M) but also on dimensionless ratios of 
kinematic variables. For example, in the case of production of heavy quarkonium at large 
transverse momentum pt, the coefficients -F„(A) depend strongly on p^/M"^. In determining 
the relative importance of the various terms in ( |6.4|) , one must take into account not only the 
size of the matrix element and the leading power of as{M) in the short-distance coefficient, 
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but also the dependence of -F'„(A) on dimensionless ratios of kinematic variables. The terms 
given explicitly in ( |6.8| ) may not be the dominant contributions to the cross sections if the 
coefficients of the matrix elements are sufficiently suppressed relative to the coefficients of 
other matrix elements that are of higher order in v. 

E. P-wave Production 

We next apply the factorization formalism to the production of P-wave quarkonium to 
leading order in f , using the lowest-lying P-wave levels of charmonium for the purpose of 
illustration. According to our factorization formula ( [6.4] ) , the inclusive production rates for 
P-wave charmonium are 

<K) = ^^^ (0|O^(^Pl)|0) + ^^ (0|08^=(^5o)|0) + 0(t;V), (6.11a) 

<xcj) = ^^ (o|of-^ep,)|o) + ^^ (oioi-e^oio) 

+ 0(t;V), J = 0,1,2. (6.11b) 



The vacuum-saturation approximation ( |6.6D can be applied to the color-singlet matrix ele- 
ments to express them in terms of vacuum-to-quarkonium matrix elements. These matrix 
elements can be expressed in terms of regularized derivatives of radial wavefunctions at the 
origin by using ( p.l9|) . Because of heavy-quark spin symmetry, they can all be replaced, 
without loss of accuracy, by their weighted average R'p. Heavy-quark spin symmetry also 



implies that the color-octet matrix elements in (|6.11| ) are proportional to 2J -|- 1, up to 
corrections of relative order f^. Thus, the P-wave charmonium production rates can all 
be expressed in terms of the two nonperturbative parameters |-Rpp and {Q\0^''{^Sq)\Q) (or, 
alternatively, the average over the P-wave states of 3/(2 J + 1) times the color-octet matrix 
elements) : 

"(h.) = ^%^ Tf/ + ^^ {0|OM'S„)|0) + 0(.=.). (6.12a) 



27rM4 



M2 



(2J + l)3Ar,Fi(3Pj) —2 ^ {2J + l)F^{^Si 

3M2 

+ 0(t;V), J = 0,1, 2. (6.12b) 



o(x.) - '■ ^ S '" ' I'^f + '--' ' motcsm 
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Note that the color-octet matrix element {0\Og''{^So)\0) in ( p.l2| ) cannot be identified with 
the decay matrix element {hdOsC So)\hc) in ( [4.7| ). 



The first application of the factorization formulas (|6.12|) was to the inclusive production 
of P-wave charmonium states in i?-meson decay p4| . The factorization formulas were given 
in the form 



K + x) 



T{b^XcJ + X) 



+ 3H',{fi)fs{b^ccCSo)+X 
Hifi[b^ccCPj)+X,fi) 

+ (2J + 1) H'sifi) fg (b -^ cci^'Si) + X 



(6.13a) 



(6.13b) 



The coefficients Fi and Fg are proportional to the production rates for on-shell QQ pairs 
in color-singlet P-wave and color-octet S-wave states, respectively. The factors Hi and Hg 
can be expressed in terms of NRQCD matrix elements divided by appropriate factors of the 
heavy-quark mass: 

1 



Hi 



(0|O^(^Pi)|0), 



3M4 



(6.14a) 
(6.14b) 



The definitions (|6.14a|) and (|6.14b|) were chosen in Ref. |2^ so that Hi and ifg would 
coincide as closely as possible with the decay matrix elements Hi and Hg. Using the vacuum- 
saturation approximation ( |6.6| ), we see that the definition of Hi given in ( |6.14a|) is equal to 
that given in ([4.9a|) , up to corrections of relative order f*^. A crude estimate for Hg{M) in 
terms of Hi is given in ( p.30|) . A similar estimate of Hg{M) in terms of Hi can be obtained 
by solving the evolution equation ( |6.5| ) and assuming that (0|C8''(^S'o; A)|0) can be neglected 
at some initial scale A = Aq. With the normalizations in ( [4.9b| ) and ( |6.14b| ), the resulting 
estimates for Hg{M) and Hg{M) are equal. However, there is no apparent rigorous relation 
between these two matrix elements. 

As we have already remarked in connection with the decay matrix elements, the factors 
of 1/M in (|6.14a|) and (|6.14b|) are more properly associated with the operator coefficients. 



since they involve short-distance physics at distance scales of order 1/A or less. Therefore, 
the factorization formulas ( |6.12| ) are preferable to the forms given in ( |6.13| ). 



In Ref. [Q, which discusses the decay of a B meson into a charmonium state, the NRQCD 
cutoff A was set equal to the scale of the large momentum transfer in the process, which 
is the bottom-quark mass ruf,. This choice of cutoff is inappropriate because the NRQCD 
evolution equation ( |6.5|) accurately reflects the behavior of full QCD only for cutoffs A that 
are less than M. That is, the NRQCD evolution equation cannot be used to sum logarithms 
of Q'^/M'^, where Q^ is the large momentum transfer in a production process. Therefore, 
a more appropriate choice of NRQCD cutoff for the process analyzed in p4| is A = rric, 



where rric = M is the charmed-quark mass. Note, however, that a change of NRQCD cutoff 
from rrih to rric does not affect the short- distance coefficients in the leading-order calculation 



presented in |24], and is, in general, insignificant numerically. 
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VII. DISCUSSION AND OUTLOOK 

The factorization approach that we have developed in this paper provides a systematic 
theoretical framework for understanding the annihilation and production of heavy quarko- 
nium. In this section, we discuss the relation between our approach and previous models 
for quarkonium production and annihilation. We also summarize the current status of the- 
oretical calculations of annihilation rates and production cross sections. 

A. Comparison with Previous Approaches 

We have presented a rigorous formalism for calculating the inclusive annihilation rates 
of heavy quarkonia. It is based on the use of NRQCD to separate the annihilation rate 
into short-distance parts, involving distance scales on the order of 1/M, and long-distance 
parts. The short- distance parts are identified with the imaginary parts of coefficients in 
the NRQCD lagrangian, and can be computed as perturbation expansions in as{M). The 
long-distance parts are expressed as matrix elements of 4-fermion operators in NRQCD and 
can be computed nonperturbatively by using lattice simulations. We have also developed 
an analogous formalism for computing inclusive production rates of heavy quarkonia in 
processes involving large momentum transfers. The cross sections are factored into short- 
distance parts, which can be computed perturbatively, and long-distance parts, which are 
expressed as NRQCD matrix elements. 

The factorization approach provides a firm theoretical foundation for calculations of the 
annihilation and production rates for heavy quarkonium. It can be used to assess the degree 
of validity and the limitations of models used in previous work on heavy quarkonium produc- 
tion and annihilation. The most thoroughly developed model for the calculation of produc- 
tion rates is the "color-singlet model" p0tpT] , p2tp3| . Most calculations of annihilation rates 
have also been carried out within this model. In the color-singlet model, the quarkonium 
state is modeled by a color-singlet QQ pair that is in the appropriate angular-momentum 
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state and has vanishing relative momentum. Nonperturbative effects are assumed to fac- 
tor into a single nonperturbative quantity that is proportional to the square of the radial 
wavef unction, or one of its derivatives, evaluated at the origin. 

The factorization formalism represents a significant advance over the color-singlet model 
in several respects. First, it provides a systematic framework for calculating perturbative 
corrections to the short-distance factors to arbitrarily high orders in Ug- The infrared di- 
vergences that are encountered at any order of the perturbation expansion can be factored 
into specific nonperturbative matrix elements. Perturbative calculations in the color-singlet 
model are based on the assumption that all infrared divergences can be factored into a single 
nonperturbative quantity. In the case of S-waves, calculations at next-to-leading order in 
as (NLO) provide empirical support for the assumption that long-distance effects can be 
factored into the quantity |-R5(0)p. Our formalism reveals that this assumption is, in fact, 
correct for any specific S-wave process in the nonrelativistic limit to all orders in as- It has of- 
ten been assumed, in addition, that the same quantity |-R5(0)p describes processes involving 
both the ^ and 1 S-wave states. Our formalism shows that this additional assumption 
is correct only up to corrections of relative order v"^. The assumption that the same quantity 
|i?5(0)p describes annihilation into light hadrons and electromagnetic annihilation also fails 
at relative order t>^, as does the assumption that the same quantity |-Rs(0)p describes both 
annihilation and production processes. In the case of P-waves, explicit calculations of the 
decay rates into light hadrons reveal that the assumption of a single long-distance factor 
|-Rp(0)p fails at leading order in as (LO) for he and Xci and at NLO for Xco and Xc2 [0- 
In the context of our formalism, these results follow simply from the existence of a second 
independent matrix element that contributes to the annihilation rates of P-wave quarkonia 
in the nonrelativistic limit. 

The factorization formalism also improves upon the color-singlet model by allowing the 
systematic calculation of relativistic corrections to annihilation and production rates. Rela- 
tivistic corrections are incorporated by including nonperturbative matrix elements that scale 
as higher powers oiv. In the case of S-waves, our formalism for computing the f ^ corrections 
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is similar at leading order in as to a model for relativistic corrections developed by Keung 
and Muzinich [Q. The major differences are that the factorization formalism provides non- 
perturbative definitions for the long-distance factors, it allows the short-distance coefficients 
to be calculated beyond leading order in a^, and it can be used to treat corrections of order 
f ^ and higher. 

Another advantage of the factorization formalism is that it provides unambiguous field- 
theoretic definitions of the long-distance factors in annihilation and production rates. This 
allows one to compute them nonperturbatively using, for example, lattice simulations of 
NRQCD. Previous approaches have relied either on determining the long-distance factors 
phenomenologically or on relating them to potential-model wavef unctions. Both of these 
approaches are of limited utility. The purely phenomenological approach can be applied 
only in situations in which the number of accurately-measured experimental observables is 
greater than the number of nonperturbative matrix elements. Potential-model estimates 
can be used for color-singlet matrix elements that have simple potential-model analogs, but 
they cannot be used for other matrix elements, such as the color-octet matrix elements 
that contribute to the annihilation of P-wave states into light hadrons at leading order in 
V. It is also difficult to gauge the accuracy of potential-model estimates in the absence 
of a rigorous connection to QCD. Since our formalism provides unambiguous definitions of 
the long-distance factors in annihilation and production processes, it allows us to quantify 
relations between these matrix elements and Coulomb-gauge wavefunctions in NRQCD. It 
also allows us to quantify the differences between matrix elements for decays into light 
hadrons and matrix elements for decays into electromagnetic final states, as well as the 
differences between annihilation matrix elements and production matrix elements. 

A final advantage of the factorization formalism is that it takes into account the complete 
Fock-space structure of the quarkonium. In the color-singlet model, the quarkonium is 
assumed to be simply a QQ pair in a color-singlet state with definite angular-momentum 
quantum numbers "^^^^Lj. However, a quarkonium also has a probability of order v"^ to be in 
a QQg Fock state, and it has probabilities of order f ^ or smaller for the higher Fock states. 
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In the case of P-waves, the factorization formahsm reveals that the QQg component can play 
just as important a role in annihilation and in production as the dominant QQ component. 
In the case of S-waves, the higher Fock states can be ignored in the nonrelativistic limit, and 
even at relative order f^, but the factorization formalism indicates that they do contribute 
at relative order v^. 

The factorization formalism for describing the annihilation of heavy quarkonia is in many 
ways similar to the operator-product-expansion formalism for calculating the inclusive decay 
rates of heavy-light mesons ||35|. These decay rates can be factored into short distance parts, 



which involve the weak decay of a heavy quark or its weak annihilation with an antiquark 
in the meson, and long-distance parts, which can be expressed as NRQCD matrix elements. 
The main difference between heavy-quarkonium annihilation and heavy-light meson decay 
is in the relative importance of the various matrix elements. Since the typical momentum of 
a heavy quark in a heavy-light meson is of order Aqcd and is independent of M, the relative 
importance of matrix elements is determined strictly by the dimension of the operator. 

Operator-product-expansion methods have also been used to treat exclusive decays of 
heavy quarkonium into light hadrons at leading order in v [Q. The NRQCD formalism 
might prove to be useful in extending such analyses to include relativistic corrections. In 
exclusive processes, a factorization theorem holds, not only for the decay rate, but also for 
the decay amplitude. Thus, just as in the case of electromagnetic annihilation, the relevant 
NRQCD matrix elements for exclusive decays are vacuum-to-quarkonium matrix elements 
of color-singlet operators of the form x^^n'4'- 

Operator-product-expansion methods have also been used in a completely different con- 
text in heavy quarkonium physics @,0. These methods have been used to treat the in- 



teractions of heavy quarkonium with light hadrons whose momenta are small compared to 



the scale Mv of quarkonium structure. Voloshin |^ has used this approach to calculate 
nonperturbative corrections to quarkonium annihilation rates that are proportional to the 
gluon condensate. In our factorization formula, the gluon-condensate contribution would 
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appear in the long-distance matrix element. In some of the cases considered by Voloshin, 
the corresponding short- distance part involves the annihilation of the QQ pair in a color- 
octet state. His approach can therefore be used as a framework for estimating the matrix 
elements of color-octet operators. 

The general factorization formula ( |6.4| ) for the production cross section of any specific 
quarkonium state H takes into account the short-distance production of color-singlet QQ 
pairs and color-octet QQ pairs in all angular- momentum states. In this respect, our ap- 
proach has some elements in common with the "color-evaporation model" for quarkonium 
production |^9|. In this model, the total inclusive cross section, summed over all quarkonium 
states H, is obtained by integrating the perturbative cross section for inclusive QQ produc- 
tion from the quark threshold 2M up to the physical threshold for the production of a pair of 
heavy-light mesons. No constraints are imposed on the color and angular momentum states 
of the QQ pair. Under the hypothesis of "semilocal duality", the nonperturbative QCD 
effects that are responsible for the formation of a color-singlet bound state containing the 
QQ pair are assumed to be negligible after one sums over all quarkonium states H. In the 
factorization approach, the nonperturbative effects are not neglected, but are factored into 
long-distance matrix elements (0|(9^|0). In the color-evaporation model, the production 
cross section for a specific quarkonium state H is obtained by multiplying the total quarko- 
nium cross section by a purely phenomeno logical fraction Jh- The relative production rates 
of different quarkonium states are, therefore, not predicted. In the factorization approach, 
the relative production rates can be calculated by using perturbative QCD, once the values 
of the dominant matrix elements (0|O^|0) have been determined. 

B. Present Status of Calculations 

The possible applications of the factorization formalism for heavy-quarkonium annihi- 
lation and production are almost limitless, since heavy quarkonia play a role in so many 
high energy processes. In order to highlight some of these applications, we discuss below 
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the present status of calculations of annihilation and production rates. 

In the case of S-wave decays, NLO perturbative corrections have been calculated for all 
the annihilation rates. In many cases, the NLO corrections are uncomfortably large. In 
order to develop a better understanding of the origin of these large corrections, it would 
be desirable to have calculations at NNLO, at least for the simplest processes ip -^ e^e~ 
and Tjc -^ '~f'~f. Relativistic corrections to the S-wave annihilation rates have been studied 



by Keung and Muzinich ||3J]. From their results, one can extract the coefficients of all the 
matrix elements of relative order v'^ at leading order in a^- A phenomenological analysis 
of the decay rates of the lowest-lying S-wave states of charmonium, including the next-to- 
leading order corrections in as{M) and the corrections of relative order f ^, is in progress 

In the case of P-wave decays, complete NLO perturbative corrections are available only 
for the electromagnetic decays Xco -^ 77 and Xc2 — ^ 77- For the decays of P-wave states 
into light hadrons, complete results are known only to order a^ |[T^. The coefficients of 
|i?pp have been calculated to order a^ 0, but they contain logarithmic infrared divergences 
that should be factored into matrix elements of the operators Os{^Sq) and Os{^Si). There 
are constants under the logarithms that should also be factored into the matrix elements. 
Unfortunately, these constants cannot be determined readily from the existing calculations. 
The relativistic corrections to P-wave annihilation rates have not yet been analyzed. 

In the case of D-wave decays, the only complete LO calculations are those for the elec- 
tromagnetic decays of the ^Di state into e"*"e~ and the ^Z^2 state into 77 0. For the decay 
of the ^D2 state into light hadrons, the coefficient of the matrix element corresponding to 
|-R^(0)p has been calculated at LO 0. For the decays of the ^Dj states into light hadrons, 
only the logarithmic infrared divergence in the coefficient of |-R^(0)p has been extracted 
Pl| . This divergence should be factored into other matrix elements that contribute to the 
annihilation rate in the nonrelativistic limit. These matrix elements can be identified by 
using the methods of Section |111 A| , and their coefficients can be calculated by using the 
methods illustrated in Appendix A. 



The status of calculations of the production of heavy quarkonium has been reviewed 



recently in Ref . |2^ , although many aspects of that review are superceded by the develop- 
ments described in the present paper. In the case of S-waves, most production processes 
have been computed only to LO. The only processes for which complete NLO calculations 
are available are if) and rj^ production in 5-meson decay [Q and inclusive rjc production 
in hadron collisions [^,^. It would be desirable to have calculations of complete NLO 
corrections for more production processes, in order to develop a better understanding of the 
size and behavior of the perturbative corrections. It is also important to calculate the rela- 
tivistic corrections, which are expected to be typically on the order of 30% for charmonium. 



Relativistic corrections have been calculated for the photoproduction of the ip [44| within 



the model of Keung and Muzinich [0]. The factorization formalism can be used to express 



those results in terms of well-defined NRQCD matrix elements. 

For the production of quarkonia at large transverse momentum px-, the contributions 
that are leading in 1/pt sometimes come from beyond leading order in the perturbation 
expansion, and they can be computed without complete calculations of the NLO or NNLO 
corrections. These contributions come from fragmentation and can be expressed in terms 
of process-independent fragmentation functions Di^niz, fJ') for a parton i with invariant 
mass /i to produce a jet containing the quarkonium H with light-cone momentum fraction 
z. The fragmentation contribution to a production cross section sometimes appears in 
a LO calculation, but it often appears first at NLO and sometimes even at NNLO. The 
fragmentation functions for producing S-wave quarkonia from the fragmentation of gluons 
| 45| and heavy quarks |^ have been calculated at LO in ag- 



For P-wave quarkonia, there are many production processes for which complete calcu- 
lations are not even available at LO. For most processes, the coefficient of |-RpP has been 
calculated |^. Complete LO calculations, including the coefficient of (0|O|^|0), are avail- 
able only for the production of P-wave charmonium in 5-meson decays |2^, T decays P7 



gluon fragmentation |^, and charm fragmentation |^D|. Relativistic corrections to P-wave 



production processes have not been studied. 
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In the case of production of D-wave quarkonia, the only perturbative calculations that 
are available are those for the coefficients of |i?^(0)p in the decay rates for Z^ —>■ Hj, where 
H represents any of the ^D2 or ^Dj states ||5^ . 



In the factorization approach, nonperturbative long-distance effects are organized system- 
atically into well-defined NRQCD matrix elements. This allows one to go beyond potential 
model estimates or phenomenological determinations of the long-distance factors. Instead, 
they can be calculated from first principles using lattice simulations of NRQCD. Such calcu- 
lations are still in their infancy. At present, the only vacuum-to-quarkonium matrix elements 
that have been calculated are (0|xVl^c), {0\x^ (^i^l'ip) , and (0|x^Dx|^c), and their analogs for 
the bottomonium system [^,^]. The only 4-fermion matrix elements that have been cal- 
culated thus far are {ric\Oi{^ So)\ric) , {hc\Oi{^ Pi)\hc) , and {hc\Os(^ So)\hc) and their analogs 
for the bottomonium system [^. Thus far, all matrix elements have been calculated only 
up to corrections of relative order f ^ and in the absence of dynamical light quarks. Pro- 
duction matrix elements are much more difficult to calculate through lattice simulations, 
unless they can be related to annihilation matrix elements through the vacuum-saturation 
approximation. 

C. Concluding Remarks 

Heavy-quark mesons have long been the best understood of hadrons. Until recently, 
our understanding has been based almost exclusively on phenomenological quark potential 
models that are motivated by QCD. Now, lattice QCD simulations are providing systematic 
analyses that are based directly upon the QCD lagrangian [^. Heavy-quark systems are 



particularly well-suited to lattice simulations, and, consequently, they are now of central im- 
portance to our exploration of nonperturbative QCD. This new role for quarkonium studies, 
as a rigorous testing ground for nonperturbative QCD, demands a much higher degree of 
rigor than was necessary in older phenomenological analyses. Approximations are necessary 
in tackling most hard problems, but it is essential in a fundamental analysis that there be 
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systematic procedures for improving the approximations. In this paper, we have developed 
a formahsm for studying annihilation decays of heavy-quark mesons that meets this stan- 
dard. With our formalism, we can improve upon the nonrelativistic quark potential model 
by including relativistic corrections in a systematic way. We can also go beyond the quark 
model to include the dynamical effects of gluons. Thus, we can, for the first time, begin 
to confront the full richness of nonperturbative QCD in analyses that are systematic and 
rigorous. 
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APPENDIX A: COEFFICIENTS OF 4-FERMION OPERATORS 



The coefficients in the lagrangian for nonrelativistic QCD can be determined by match- 



ing scattering amphtudes in NRQCD with those in full QCD \T^\. In this Appendix, we 
use these techniques to determine the coefficients for some of the 4-fermion operators that 
contribute to quarkonium annihilation rates. In Section |A1|, we illustrate the method by 



calculating the coefficients of dimension-6 operators to order as- In Section |A2|, we apply the 



method to calculate the imaginary parts of the coefficients of dimension-6 and dimension-8 
operators to order a^. In Section [Aij|, we demonstrate how the imaginary parts of some of 



the coefficients can be extracted at next-to-leading order from existing calculations of the 
decay rates of bound states. We also record coefficients that can be extracted from existing 



calculations in the literature. Finally, in Section |A4| , we give the corresponding coefficients 
for electromagnetic annihilation rates. 

1. Coefficients at Order Ug 

We wish to determine the coefficients of the dimension-6 and dimension-8 4-fermion 
operators at order a^ by using the matching condition ( |2.18| ). We consider QQ scattering 
amplitudes, with the momenta of the heavy quarks and antiquarks small compared to the 
heavy quark mass M. In full QCD, there are two Feynman diagrams for QQ scattering at tree 
level. The gluon exchange diagram in Fig. ^(a) is also present in NRQCD. The annihilation 
diagram in Fig. Kb) is not present in NRQCD, so its effects must be reproduced by adding 
4-fermion terms to the effective lagrangian. We calculate the annihilation contribution to 
the amplitude for QQ scattering in the center of momentum frame. We take the incoming 
Q and Q to have momenta p and — p, while the outgoing Q and Q have momenta p' and 
— p'. By conservation of energy, we have |p'| = |p| = p. 

The scattering amplitude (T-matrix element) in full QCD corresponding to the diagram 
in Fig. Kb) is 
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M 



6(b) 



TTOs 



u{p')rT''vi-p') v{-ph,T^uip), 



(Al) 



where E = \/M'^ + p^. We have suppressed the color indices on the 4-coniponent Dirac 



spinors. Following Ref. |T^, we express the 4-coniponent Dirac spinors in the Dirac repre- 
sentation in terms of 2-coniponent Pauli spinors via the substitutions 



u[p) 



v[-p) 



lE + M 



( 



2E 



e 



\ 






lE + M 



2E 



\ E+M>' J 
E+M I 

V ^ ; 



(A2a) 



(A2b) 



where ^ and t] are 2-component spinors with suppressed color indices. The Dirac spinors 
u{p') and f (— p') have similar expressions in terms of Pauli spinors ^' and t]'. The spinors 
( |A2a|) and ( |A2b|) represent fermion states with the standard nonrelativistic normalization. 
Expanding to second order in the velocity v = p/E, we find that the annihilation contribution 
to the scattering amplitude ( |A1| ) from full QCD reduces to 



M 



6(b) 






{1-v 



2\t'\ 



)E,'^ (tT^t]' ■ r]^(TT^^ {v'v^ + v%'^)^'^a'T''r]' 7]'^a^TX , (A3) 



where v = p/E and v' = p' / E. It is convenient to suppress the spinors and write the above 
matrix element as a direct product of color matrices multiplied by a direct product of spin 
matrices: 

1 



M 



6(b) 






(ja ^ rj,a^ 



'l~v^)a' ®a' 



-{y'v^ + v'%'^)a' ® a^ 



(A4) 



One can read off the dimension-6 term in the scattering amplitude in terms of the pa- 
rameters of NRQCD by substituting ^, ^' , r]' , and t]^ for ip, ip^ x, and x^ iii the effective 
lagrangian (p.9|): 



Md=, = ]^ (1 ® 1) [/iC-^o) 1 ® 1 + /i('5i) (y' ® o' 

+ ^{T^® T«) [fsCSo) 1 ® 1 + fsi'Si) a^ ® 



a 



(A5) 



Comparing ( [A4D and ( |A5| ) , we find that only one of the four terms in ( |2.9D has a nonvanishing 
coefficient at order a,: 
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fsCSi] 



-TvasiM). 



(A6) 



The coefficients of the remaining three terms in ( |2.9| ) are of order a^. 

To determine the dimension-8 coefficients, we need the scattering amphtudes from the 
term ( p.ll| ) in the effective lagrangian: 

1 



M, 



M2 



;i®i) 



/lOPOv'-vl®! + ^^£^^^^1±A^ v' • V a^ ® a^ 



f,{^P,)-f,{^P, 



-\ V ■ cr (g) V • cr 



+ gi[ Sq)v 1 O 1 H V 



V • cr (g) V ■ <T 



o-*(g)cr' 



+ 



^iC^l,'/^! 



[v^v^ + f '*t;'-'')cr* ® 0"^ 



+ .... 



(A7) 



There are similar terms with color structure T" ®T'^ and coefficients /s and ^fg- Comparing 
with (|A4|) , we find that 



An 
gsCSi) = —as{M) 

gsCSu^Di) = nas{M). 



(A8a) 
(A8b) 



The color-singlet coefficients gi and the remaining color-octet coefficients vanish at this order 
in a,(M). 

2. Imaginary Parts at Order a^ 



We now turn to the calculation of the imaginary parts of the coefficients at order a^. 
They can be determined by matching the imaginary parts of QQ scattering amplitudes in full 
QCD and NRQCD in accordance with ( p.l8|) . In full QCD, the annihilation contributions 
to the imaginary parts at order a^ come from the one-loop diagrams in Fig. |^. We will 
determine the imaginary parts of the coefficients of the dimension-6 operators in ( p.9| ). We 
will also determine the imaginary parts of the coefficients of the dimension-8 operators that 
contribute to the annihilation of P-wave states at leading order in v and S-wave states 



through relative order v'^. The dimension-8 terms in the lagrangian are given in ( p.lll) . To 
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determine the coefficients, we consider QQ scattering in the center of momentum frame, 
with the Q and Q momenta small compared to the heavy-quark mass M. We calculate the 
imaginary parts of the diagrams in Fig. |^ in Feynman gauge. After making the substitutions 
( |A^ ) for the Dirac spinors, we expand in powers of the velocity v. 

Below, we list the results for the imaginary parts of each of the diagrams, suppressing 
the spinors as in ( [A4| ) . The diagrams in Fig. |^(a) and |^(b) yield, on expansion through 
second order in the velocity f , 



ImA^7ra^ = ?^ {T^T^ ® T^T") 



'7(a) 



2M2 



4 1 

l--v^ + -^r■V] 1 ® 1 



V3 5 5 J 

+ (-v'v'^ + —v''v^ - —{v'v^ + v''v'^) ) a* O (T^ 

\0 -LO 0\J 



ImM-i - 



17(b) 



2M2 



^rparpb ^ rparpb^ 



l--v^ - -^r -V] 1 ® 1 



1 1 2 2 

V V ■ V I cr" (X> cr" 

3 5 5 

V5 15 30 



(A9a) 



(A9b) 



The color matrices can be simplified as follows: 



rparph ^TA n~^brjia 



rparpb ^ rparph 



r^ j\t2 2 

^ 1®1 + ^T^T"®T", 



2N, 
Cf_ 

2Nr 
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2N, 

1 



Nr 



T^^T 



(AlOa) 
(AlOb) 



where Cp = (^c -'-)/(2^c) is the Casimir for the fundamental representation. The diagrams 
in Fig. 0(c) and |^(d), which involve the triple-gluon vertex, yield 



ImM 



7(c) 



NcTTal 

"qap 



{T" ® T" 



10 



1 1 

5 2 ' 



(Alia) 



TV TTcy'^ 



1 - —v^ ] a' ®a 



-v'v^ + -v''v'^ ] a' (g) a^ 



(Allb) 
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From the gluon-loop diagram in Fig. |^(e) combined with the associated ghost loop diagram 
in Fig. 0(f), we obtain 

1 



Im (M7{e) +M 



7(f) 



12M2 ^ ' 



:i - v^W ® a' 



-(t;V +v v^)(7'- (g)o- 



(A12) 



The quark loop diagram in Fig. ^(g) gives 



ImM 



7(g) 



"gap 



(T" ® T") 



{l-v^)a'®a' {vV +v"v'^)a'®a^ 



(A13) 



Adding the amplitudes ( |A9|) , and (|A11|) - ([A13| ) and comparing with (|A5|) , we can read off 
the imaginary parts of the coefficients of the dimension-6 operators: 



ImhCSo) 






Im/8(%) = ^^^^^V^«^(M), 



lmfsCS^ 



Tcn 



4Ar, 

/ ^,2 



6 



ai{M). 



(Al4a) 

(Al4b) 
(Al4c) 



The imaginary part of the coefficient /i(^5'i) vanishes at order a"^. Comparing with 
we see that the coefficients of the color-singlet dimension-8 operators are 



Im/i(3Po) 
Im/i(3P2) 
Im 5-1 (^5*0) 



SvrC 



F „2 



2Nr 



af(M), 



27tC 



F „2 



SAT, 



«:(M), 



27tCp 9 / , ^n 



(Al5a) 
(Af5b) 
(A15c) 



The imaginary parts of the coefficients /i(^Pi), /i(^Pi), gi{'^Si) and (71(^51, '^Di) vanish at 



order a^. 



3. Imaginary Parts at Higher Order in as 



According to the matching condition ( |2.18| ), the coefficients of the 4-fermion operators 
can be computed at next-to-leading order in as by calculating scattering amplitudes at next- 
to-leading order in full QCD and equating them to the scattering amplitudes in NRQCD, 
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calculated to next-to-leading order in as- For some of the 4-ferniion operators, the imaginary 
parts of the coefficients can be extracted at next-to-leading order in a^ from calculations of 
heavy- quarkonium annihilation rates that already exist in the literature. 

As an illustration of our factorization approach, we discuss in detail the calculation of 
Im/i(^S'o) at next-to-leading order in ctg. In order to determine Im/i(^S'o), we consider the 
matrix element M for the forward scattering of a QQ pair above threshold in a color-singlet 
spin-singlet state with relative velocity 2v. The imaginary part of M. can be expressed as 
a sum over cuts through the Feynman diagrams for forward scattering. The annihilation 
contribution to IniA^ is the sum over cuts through gluon and light quark lines only. It has 
been calculated in full QCD at next-to-leading order in Ug and in the limit f — ^ by Barbieri 
et al. |E|: 



, ,, TTCFal{2M)\^ r/vr^ vr^ \^ /l99 IStt^X 



r' 



TT 



(A16) 



where Ca = N^. is the Casimir for the adjoint representation, as{M) is the QCD coupling 
constant in the modified minimal subtraction (MS) renormalization scheme for QCD with 
Uf fiavors of light quarks, and M is the perturbative pole mass of the heavy quark. We have 
corrected apparent errors in Ref. |^ of 2/3 in the overall coefficient and 1/2 in the coefficient 
of the TT^/f term. The next-to- leading order correction contains a Coulomb singularity 
proportional to 1/f , which gives an infrared divergence in the limit f —* 0. 

In order to determine Im/i(^S'o), we must calculate the corresponding contribution to 
ImA^ in NRQCD at next-to- leading order in a^. The relevant Feynman diagrams are 
shown in Fig. p. They contain a 4-fermion vertex that corresponds to the term ip'^xx^'4' iii 
the effective lagrangian. The annihilation contribution is the sum over all cuts that pass 
through that vertex. The Cutkosky rules specify that a cut passing through the 4-fermion 
vertex is computed by taking the imaginary part of the coefficient /i(^5'o) and complex- 
conjugating the part of the diagram to the right of the cut. The incoming and outgoing 
states consist of a QQ pair in a color-singlet spin-singlet state with relative velocity v ^ 0. 
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Explicitly, the in state is 

\QQ) = Y,^5'^ Yl -^{^cT2)mm'\Q{p,m,^)Q{-p,m',J)) , (A17) 

where p = Mv is the momentum of the quark in the center of momentum frame. The quark 
states \Q{p,m,i)) with momentum p, spin quantum number m = ±1/2, and color i have 
the standard nonrelativistic normalization: 

{Qip',m',j)\Q{p,m,z)) = 6'^^"^'{27c)'6\p' - p) . (A18) 

For the leading-order diagram in Fig. ^(a), the cut through the 4-fermion vertex gives simply 

the imaginary part of the coefficient /i(^5'o)/M^: 

2N,lmM'So) 
ImA48(a) = j^^ • (A19) 

It is convenient to calculate the next-to-leading order diagrams in Figs. H(b) and H(c) by 
using Coulomb gauge, since then only Coulomb exchange contributes in the limit v —* 0. 
For the diagram in Fig. H(b), we obtain 

8(b) 






1 1 

(A20) 



^ + go - (P + q) V2M + ie E - go - (p + q) V2M + ie ' 
where E = p^/2M. After using contour integration to integrate over the energy go of 
the exchanged gluon, we find that the contribution reduces to an integral over the gluon's 
3-momentum: 

V /8(b) ^ ^ (27r)^q^q^ + 2p-q-«e 

The integral is infrared divergent, and can be regularized by using dimensional regularization. 
The integral over q is analytically continued to D = 3 — 2eiR spatial dimensions. Evaluating 
the regularized integral in ( |A21 ), we obtain 



j^_^, _ 2iV,Im/iC^o) nC^a, 



8(b) 



M2 Av 



i I , , , , 2Mv\ 
1 — + log(47r) - 7 - 2 log 

(A22) 
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where 7 is Euler's constant and /iiR is the arbitrary regularization scale introduced with 
dimensional regularization. The logarithmic infrared divergence appears as a pole in ejR. 
The subscripts IR on e and fi serve as a reminder that they are associated with infrared 
divergences. Note that ( |A22| ) is complex valued. The imaginary part of ( |A22| ) arises because 
it is possible for the incoming quark and antiquark to scatter on-shell before annihilating at 
the 4-fermion vertex. After summing over all diagrams, one must, of course, obtain a real 
result for Im^Vl. The diagram in Fig. |Kc) is evaluated in the same way as Fig. ||(b), except 
that the Cutkosky cutting rules require the complex-conjugation of the part of the diagram 
that involves the Coulomb-gluon exchange. The result is 



ImM 



2iV,Im/i(i5o) nCpas 



8(c) 



M2 



4v 



1 + - — + log(47r) - 7 - 2 log 



(A23) 



Note that the imaginary part of ( |X23|) cancels that of ([A22|) . Adding ([Al9|) , ([A22D , and 
( |A23|) , we obtain the complete result for ImA^ through next-to-leading order in a^: 



ImM 



2N,lmM^So) 
M2 



1 + ^^c,^ 

2v IT 



(A24) 



Comparing ([A16|) and (|A24|) , we can read off the imaginary part of /i(^S'o) through 
next-to-leading order in a^: 



ImfiCSo 



nC 



P ^.2 



2N,. 



<(2M)< 



{-I(t-) 



c. 



199 IStt^ 



18 24 



C, 



9 



^/ 



TT 



(A25) 



Note that the factorization approach reproduces the standard prescription of simply drop- 



ping the 1/v terms in the perturbatively calculated annihilation rate |5^. The factorization 
approach puts this prescription on a rigorous footing, and makes it clear how to extend the 
calculation systematically to higher orders in a^ and in v. 

In ( [A25| ), as{M) is the MS coupling constant with renormalization scale M. If we make a 
different choice for the renormalization scale /i of ots{^J')■, then we must differentiate between 
the MS coupling constant aT^ (/i) for full QCD with nj flavors of light quarks and a heavy 
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quark and the corresponding coupling constant aT^\^) for only n/ flavors of light quarks, 
which is the appropriate running coupling constant below the heavy-quark threshold. These 
coupling constants satisfy the matching condition [|^ a^* (M) = a^^ (M) + 0(0;^). If 
we wish to use a different renormalization scale /i 7^ M for a^ in (|A25 ), then we must make 
one of the following substitutions: 






1 + /5olog^^ + 0{al) 

M TT 



^ + (*4)'°«^v + «("•' 



(A26a) 
(A26b) 



where /3o = (33 — 2nf)/6 is the flrst coefficient in the beta function for QCD with Uf flavors 
of light quarks: /i((i/(i/i)as(/i) = —/3oal/n + . . .. 

The coefficient of the operator Oii}So) in the NRQCD lagrangian is fi(}So)/M^, and 
the perturbation series for /i(^S'o) depends on the deflnition of the heavy quark mass M. 
The order-a^ correction in ( |A25| ) corresponds to the choice M = Mpoie, where Mpoie is the 
perturbative pole mass, i.e., the location of the pole in the heavy-quark propagator in per- 
turbation theory. An alternative choice is the running mass M{fi) in the MS renormalization 

m 



scheme. Its relation to the pole mass through order as is 



l+|l + 5,ogA)c,^+0(„2) 



(A27) 



Throughout this paper, we will adopt the choice M = Mpoie for the heavy quark mass in 
the coefficient fn/M'^"~'^ of a 4-fermion operator with naive scaling dimension dn- 

We can obtain the imaginary part of the coefficient fi{^Si) through next-to- leading order 
in as from a calculation by MacKenzie and Lepage of the annihilation decay rate of the J/i/j 
or T 1^. Their published result is given explicitly only for Nc = 3, but one can insert the 
appropriate color factors in the various classes of diagrams and obtain the result 



imhCs,] 



(7r^-9)(Ar^^-4)C^ 

54Nr. 



al{M) 



+ -Q'{T.Qi) 



a 



1 + (-9.46(2)Cp + 4.13(17)C^ - 1.161(2)r2^ 



13 as 
-r'^F — 

4 TT 



a. 



vr 



(A28) 
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where Q is the electric charge of the heavy quark {Q = +2/3 for the charmed quark and 
Q = —1/3 for the bottom quark) and the Qi, i = 1, . . . ,nf, are the electric charges of the 
light quarks. The perturbative correction in the first term on the right side of ( [A28D was 



calculated by Mackenzie and Lepage |^. The term proportional to a^ is due to annihilation 
of the QQ pair into a virtual photon, which then decays into light hadrons. The order-a^ 
correction can be calculated as the sum of two terms: —ACpas/TT, which is the order-a^ 
correction to the rate for ip -^ e~^e~, and 3Ci?as/(47r), which is the order-a^ correction to 
the rate for 7* -^ qq. For completeness, we also give the coefficient analogous to (|A28|) for 
the decay of the ip into a photon plus light hadrons: 



2 _ Q\r^n'^n, _ r n, ^ 

(A29) 



2('^-^)^-«!^a2(M) 



QJ 

l + (-9,46(2)C, + 2,75(ll)C.-0,774(l),v)^ 



3N, 
Calculations of the annihilation rates of P-wave states were carried out through order a^ 



s 



by Barbieri and collaborators [[7|,^,PT|. They calculated only the coefficients of |-Rpp in ([4.71). 



These coefficients contain logarithmic infrared divergences that should be factored into the 
color-octet matrix elements, along with associated constants that can be determined from 
calculations in NRQCD. In Ref. [[7|, the logarithmic infrared divergences in Im/i(^Po) and 
Im/i(^P2) were cut off by taking the heavy quark and antiquark off their mass-shells and 
below threshold, in which case the infrared divergence manifests itself as a logarithm of the 
binding energy. In order to extract NRQCD coefficients, it might be necessary to repeat 
the next-to-leading order calculations in Ref. 0] using on-shell scattering amplitudes and 
dimensional regularization of the infrared divergences in order to maintain gauge invariance. 

4. Coefficients of Electromagnetic Operators 

The calculation in Section |A2| can be easily modified to give the imaginary parts of the 



coefficients of the electromagnetic 4-fermion operators at order a^ and at leading order in 
as- The Feynman diagrams in Figs, ^(a) and ^b) yield imaginary parts that correspond 
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to annihilation into two photons. These imaginary parts can be obtained from (|A9|) by 
replacing the color matrices T" by the unit color matrix and by substituting a^ — > Q^a, 
where Q is the electric charge of the heavy quark: Q = +2/3 for the charmed quark, and 
Q = —1/3 for the bottom quark. The sum of the 2 diagrams yields 



4^,2 



Im iMg(^)+M 



'9(b) 



nQ a 



:i^v. 



1 v^] 1 (g) 1 H V ■ V 0-* ® cr 



'2 11 • 

+ ( -v'v'^ + —v"vJ ) a' ® oJ 
.5 15 



(A30) 



Comparing to the NRQCD scattering amplitudes analogous to ([A^), we find that the only 
nonzero coefficient for the dimension-6 operators is 



Im/^^("^S'o) = -kQ^c?. 



(A31) 



Comparing to the NRQCD scattering amplitudes analogous to (|A7|) , we can read off the 
nonzero coefficients of the dimension-8 operators: 



Im/^^(3Po) = 37rQV, 
Im/^^(^P2) = —^, : 



Im 5-^,^,(^50) 



4^,2 



47rQ4a 



(A32a) 
(A32b) 

(A32c) 



The diagram in Fig. ^(c) yields an imaginary part that corresponds to the annihilation 
into lepton pairs. The imaginary part can be obtained from ( |A13| ) by replacing T" by the 
unit color matrix, and by substituting {nf/2)as — > —Qa. The resulting matrix element is 



2^,2 



ImA^ 



9(c) 



ttQ a 
3M2 



:i®ii 



:i - v^)a' a' 



;(t;V +v''v'^)a' 0a^ 



(A33) 



Comparing to the NRQCD scattering amplitudes analogous to 
nonzero coefficient of the dimension-6 operators is 

Im/ee( Si) = . 



we find that the only 



(A34) 



Comparing to the NRQCD scattering amplitudes analogous to (|A7D, we can read off the 
nonzero coefficients of the dimension-8 operators: 
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Im5'ee(^5'i^ 



2^,2 



9 



(A35a) 
(A35b) 



Several of the electromagnetic coefficients can be determined through next-to-leading 
order in a^ from calculations that are available in the literature. The annihilation rates for 
Vc, Xco, and Xc2 into two photons have been calculated through next-to-leading order in a^ 
by Barbieri et al. [^0. The corresponding coefficients are 

^2 \ 



Im/^^(^S'o) ='kQ 



W 



^Mt-K^v 



Im/^^(3Po) =37rQV 



,4 3J ^TT 



1 - 4C, 



a.o 



vr 



(A36a) 
(A36b) 
(A36c) 



The rate for ip -^ e^e is known through next-to-leading order in as [^Z 

7rQ2a2 



Im/eer^l) 



4Cf— 



TT 



(A37) 



Finally, Mackenzie and Lepage |0 have calculated the rate for if) — > 777 to next-to-leading 
order in as- The corresponding coefficient is 

4(7r2 - 9)Q6a3 " 



Im/3,('5i; 



9 



9.46(2)C^ 



as 



TT 



(A38) 
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APPENDIX B: EVOLUTION OF 4-FERMION OPERATORS 



As we mentioned in Section |I11E|, loop corrections to the 4-ferniion operators in the 



NRQCD lagrangian are, in general, ultraviolet divergent, and, therefore, must be regularized. 
One can remove power divergences, either by employing a mass-independent regularization 
scheme, such as dimensional regularization, or by making explicit subtractions. Once this 
has been done, the 4-fermion operators satisfy simple evolution equations of the form ( ^.22p . 
The evolution equation for an operator On with naive scaling dimension dn involves only 
operators Ok with dimensions d^ > dn- The coefficients jnk in the evolution equation can 
be computed as power series in a^. For dk = dn, the coefficients 'jnk are at most of order 
a1, because the logarithmic ultraviolet divergences at order Og come only from corrections 
of relative order v^, which correspond to operators Ok of dimension (i„ + 2 or larger. In this 
Appendix, we compute at order a^ the coefficients of the dimension-8 operators that appear 
in the evolution of the dimension-6 4-fermion operators. 

1. Heavy Quark Self-energy 

In order to illustrate the methods that are used to calculate the coefficients in the evo- 
lution equations, we first calculate the self-energy of the heavy quark in NRQCD through 
order a^. From this calculation, we determine the relation between the perturbative pole 
mass Mpoie and the mass parameter M in the NRQCD lagrangian, and we extract the residue 
Z{p) of the pole in the heavy-quark propagator through order a^f ^. The residue is given by 

Z{p)-' = 1 - — (i? = pV2Mpoie,p), (Bl) 

where T,{E,p) is the self-energy correction. To determine Z{p) to order a^ and to order v"^, 
we must calculate the self-energy correction E that arises from the one-loop diagrams in 
Fig. |iy. We calculate these diagrams in Coulomb gauge, because it facilitates the extraction 
of the dependence on v. The seagull diagram in Fig. ^(b) gives only power ultraviolet 
divergences, which are subtracted as part of the regularization scheme. Interactions from 
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'^bilinear also need not be included, because the terms of order agv"^ that they produce are all 
proportional to power ultraviolet divergences. 

The contribution to the self-energy from the diagram in Fig. 0(a) is 

/d^q 1 / 1 p^ — fp • q)^ \ 

12^0^ i^-go-(p-q)V2M + 2el,^ + M^ql - ci^ + it))- 

(B2) 

The integral of the term containing 1/q^, which comes from Coulomb exchange, gives rise 
to an ill-defined power divergence, which can be dropped. After using contour integration 
to integrate over the energy go of the gluon, we find that the contribution reduces to 

^ 2'nCFas f d^q 1 P^ - (p • q)^ .^.on 

^ '^' M2 J (27r)3 q E-q-{p- q) V2M + le ' ^ ^ 



In order to identify the power divergences in ( |B3| ), we expand the denominator in a Taylor 
series in 1/M: 



2nCFas f <Pq p^ - (p ■ g) ^ ^ , i^ - pV2M + (2p ■ q - g2)/2M , 
^^^'^) = -^1^7(2^ q- i' + -q + 



(B4) 



Setting E = p'^/2M, we find that every remaining term in the integrand in (|B^ ) yields 
a power divergence, which is subtracted in our regularization scheme. Thus, after regu- 
larization, the self-energy vanishes on the energy shell, and there is no correction to the 
energy-momentum relation E = p'^/2M. In full QCD, the energy-momentum relation de- 
fined by the pole in the perturbative heavy-quark propagator is E'^ = p^ + ^poie- Matching 
the coefficients of p^ in these energy-momentum relations, we obtain 



M = Mpoie (1 + 0{ai)) . (B5) 

Thus, through order a^, the mass parameter M in the lagrangian for NRQCD can be 
identified with the perturbative pole mass. 

We proceed to compute the residue of the pole in the heavy-quark propagator, which is 
given by (|B 1|) . After we subtract the power divergences, the only term remaining in ( p4D 
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that contributes at order UgV^ is the term proportional to E — p^ /2M. Consequently, the 
expression for the residue is 

Imposing a momentum cutoff |q| < A on the magnitude of the gluon momentum and keeping 
only the logarithmic ultraviolet divergence at order a^, we find that the residue Z{p) is 

Z{v) ~ 1 ^ V , (B7) 

where v'^ = p'^/M'^. 



2. One-loop Ultraviolet Divergences 

The coefficients in the evolution equation for the operator Os(^Sq) = ip'^T'^xX^T^ip can 
be determined at order Ug by computing one-loop corrections to scattering amplitudes in 
NRQCD that involve this operator. We consider the amplitude for the scattering of a QQ 
pair with momenta p and — p into a QQ pair with momenta p' and — p'. We use the 
compact notation with suppressed Pauli spinors that was introduced in Eq. ( |A4| ) . The 



matrix element corresponding to the leading-order diagram in Fig. ^(a) is then written 

M = {T^^T") (1®1) , (B8) 

where the first factor gives the color structure and the second factor gives the spin structure. 
The one-loop correction to the matrix element is given by the sum of the contributions of the 
10 diagrams in Figs. |lT](a)-|ll|(j). We wish to calculate the terms in these contributions that 
are proportional to log A, where A is an ultraviolet cutoff. For higher-order calculations, it 
might be wise to impose the cutoff by using dimensional regularization, in order to maintain 
gauge invariance, but for our purposes it is sufficient to impose a cutoff on the magnitude 
of the gluon 3-momentum: |q| > A. 

The 4 diagrams in Fig. ^(a)-|TT](d) are self-energy corrections to the external quark 
lines. Each diagram contributes \fZ — l times the leading order amplitude in (|B^ ). Using the 
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expression ( p7|) for the renormalization constant Z, we find that the sum of the contributions 
of the 4 diagrams is 

-Mn(a-d) - -l^^il^ v' iT'^ ^ T'^) (1^1). (B9) 

The diagram in Fig. |TT](e) represents the exchange of a transverse gluon between the incoming 
quark and antiquark. (The exchange of a Coulomb gluon does not lead to an ultraviolet 
divergence.) This diagram yields the contribution 

.•47ra. (rj^a ^ rrbrrarrb\ ^T ^ 1^ f d*<l P^ - (P ■ q)^ 

(27r)4 q^-cf + ie 

^ (BIO) 



A^n(e) = .^(T«®T^T"T^) il^l) J^ 



E + qo-{p + q) V2M + te E ~ go - {p + q) V2M + ie ' 

where E = p^ /2M. We integrate over the energy go of the exchanged gluon and identify 
the power divergences by expanding the denominators in a Taylor series in 1/M. Keeping 
only the term that gives a logarithmic ultraviolet divergence, we find that the contribution 
reduces to 



M, 



The integral is the same as in ( |B6D . The diagram in Fig. |TI](f) gives an identical contribution: 

-Mll(e) ^ -Mll(f) ^ J_^I^^2 ^j.a^j.tj.aj.1,^^ (^ ^ ^^ _ (3^2) 

The diagrams in Figs. pJ](g)-[Tll(.i) involve the exchange of a transverse gluon between 
initial and final quark or antiquark lines. (The exchange of a Coulomb gluon leads to a 
vanishing contribution.) These diagrams are evaluated in the same way as those in Fig. |ll](e). 
The results are 

-Mii(g) ^ -Mii(h) ^ ^^^^vv' {t^T'^T'^T") (1®1), (B13a) 

-Mii(i) ^ -Mii(j) ^ ^^^^vv' (T'^T'^T'T'^) (1®1). (B13b) 



107 



The color factors in ( pi3|) can be simplified by using the identities in ( [A10|) . Adding up the 



results for the diagrams in ( [B9| ) and ( |B12D -( [Bl3| ), we find that the sum of the logarithmically 
divergent terms of order a^f ^ is 



+ 



{N^ - 4) V ■ v' - {N^ - 2) v"^] (T" ® T") (1 ® 1) . (B14) 



The logarithmically divergent part of the diagrams for scattering through the color-singlet 
operator Oi{^So) = 'i/'^xxV can be obtained from the expressions ( p9D and (pi2| )- (|B13b|) 



simply by replacing the color matrix T"- by the unit matrix 1. Adding up these contributions, 
we obtain 

MiCSo) ~ ^^ii^ (v ■ v' (T" ® T") - Cpv^ (1 ® 1)) (1 ® 1) . (B15) 

The ultraviolet divergent parts of the matrix elements Aigi^Si) and Aii{^Si), which corre- 
spond to scattering through the spin-triplet operators Os{^Si) and Oi(^S'i), can be obtained 
by replacing the spin factor 1 (g) 1 by a* ® o"* in ([B14| ) and ( [B15| ), respectively. 

3. Evolution Equations 

The logarithmically divergent contributions to the scattering amplitudes in Section P2 



can be expressed to leading order in v as the matrix elements of dimension-8 operators. 
Differentiating the operator equation corresponding to ( |B15| ) with respect to A, we obtain 
the evolution equation for the operator Ci(^S'o): 

By differentiating the operator equation corresponding to (pi4|) with respect to A, we obtain 
the evolution equation for the operator Os(}So): 

A— 08( So) - ^^^^^-^, 0,{ P,) + 3^^^^, Os{ Pi) - 3^^^^2 ^8( So) . 

(BIT) 
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The evolution equations for the corresponding spin-triplet operators can be obtained from 
( [B17D and ( |B16| ) simply by inserting a* between ip"^ and x ^i-nd also between x'^ ^-iid tp. It 
is convenient to express the resulting operators in terms of the combinations that appear in 
( p.ll| ) by using the identity 



D'a^ ® D'a^ = -D • cr ® D • cr + -(D x cr) 



The resulting evolution equations are 



(D X (tY + D^'a^^ ® D^'a^^ . (B18) 






5«., 



37rM2 
ACpas 

37riV^M2 



SCpCis 



(OsCPo) + OsCPi) + OsCP2)) - ^^^PiC^i, 



37rM2 



(B19a) 



Oii'Po) + Oii'Pi) + OiCP2: 



+ 



2{N^ - 4)a. 
SirNrM^ 



(OsCPo) + OsCPi) + OsCP2 



37rNrM^ 



VsCS^). (B19b) 



The evolution equations for electromagnetic operators can be calculated in the same way, 
except that there are no contributions from diagrams such as those in Figs. pl](g)-|TTl(j), which 
involve exchange of gluons between initial and final quark lines. The evolution equations 
for the dimension-6 electromagnetic operators can be obtained from (pi6|) and (pi9a|) by 
dropping the color-octet terms on the right sides and inserting vacuum projections: 



A 



d 



(^^x|o)(o|xV 



A d 



^to-xio) • (0|xV 



ACpas 
37rM2 
ACpois 



\p^x\0){0\x\-fif^ + h.c. 



V^W|0)-(0|xV(-iD)2V; + h.c. 



(B20a) 
(B20b) 
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FIGURES 
FIG. 1. Example of a diagram that contributes to the quarkonium annihilation rate at order 

a^. The three cuts of the diagram participate in a KLN cancellation. 

FIG. 2. Schematic representation of the topological factorization of the rate for quarkonium 
annihilation. The short distance part is represented by the circle labelled H. The quarkonium 
wavefunctions are represented by the shaded ovals. The wavefunctions can be connected by light 
partons, such as the two gluons that are shown explicitly. Soft gluon interactions between the light 
partons are represented by the circle labelled S. 

FIG. 3. Example of a Feynman diagram for quarkonium annihilation at order a^. The shaded 
ovals represent the quarkonium wavefunctions. 

FIG. 4. Examples of real-gluon emission in quarkonium decay at order a^. The shaded ovals 
represent the quarkonium wavefunctions. 

FIG. 5. Examples of virtual- gluon emission in quarkonium decay at order a^. The shaded ovals 
represent the quarkonium wavefunctions. 

FIG. 6. Feynman diagrams for QQ scattering at leading order in q^. 

FIG. 7. Feynman diagrams that contribute to the imaginary part of the amplitude for QQ 
scattering at order a^. 

FIG. 8. Feynman diagrams in NRQCD for the scattering of a QQ pair in a color-singlet ^^o 
state through the operator ip^xx^i^- 

FIG. 9. Feynman diagrams that contribute to the imaginary part of the amplitude for electro- 
magnetic QQ scattering at order a^. 

FIG. 10. Feynman diagrams in NRQCD for the self-energy of a heavy quark at order ag- 
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FIG. 11. Feynman diagrams in NRQCD that contribute to the evolution of an S-wave 4-ferniion 
operator, such as ^''"T^XX^^"^' or V'^XX^V'- 
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TABLES 



Operator 



Estimate 



Description 



Ois 


V 


^ 


{Mvfl^ 


X 


{Mvfl^ 


Dt (acting on V' or x) 


Mv^ 


D (acting on 'ijj or x) 


Mv 


9^ 


M2«3 


ffB 


M^v^ 


g(j) (in Coulomb gauge) 


Mv^ 


gA (in Coulomb gauge) 


Mv^ 



effective quark-gluon coupling constant 
heavy-quark (annihilation) field 
heavy-antiquark (creation) field 
gauge-covariant time derivative 
gauge-covariant spatial derivative 
chromoelectric field 
chromomagnetic field 
scalar potential 
vector potential 



TABLE L Estimates of the magnitudes of NRQCD operators for matrix elements between 
heavy-quarkonium states in terms of the heavy-quark mass M and the typical heavy-quark velocity 
V. The estimates shown apply to matrix elements in a quarkonium state \H) whose position is 
localized to a region of size 1/Mv or less. If the states are normalized to {H\H) = 1, then the 
product of the magnitudes of the operators gives the magnitude of the matrix element. (In order 
to obtain estimates for matrix elements between momentum eigenstates that are normalized to 
{H\H) = V, where V is the volume of space, one should multiply the estimates for localized states 
of unit norm by (Mv)^^.) 
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, In tills erratum, we clarify the velocity-scaling rules for those NRQCD matrix elements 

h>. ■ 

Q\^ I whose leading contributions come from \QQg) Fock states that can be reached through a 

cn ' 

en ' spin-flip transition from the dominant Fock state. A correct accounting of these spin-flip Fock 

o 

0^ 



states leads to revisions of the error estimates In several equations In the paper. In addition, 
we emphasize that the velocity-scaling rules should be used to estimate the probabilities of 
higher Fock states, rather than their amplitudes. We also correct some typographical errors. 



Oh! 
^ ' • Throughout the paper, phrases of the type "amplitude of order f "" should be replaced 



with "probability of order t>^"". The reason Is that the probability of a \QQg) Fock 



c^ ' state Is the square of the amplitude Integrated over the phase space of the particles. 

Some of the dependence on v arises from the Integration over the phase space of the 
gluon. 

Throughout the paper, one should keep In mind that the velocity expansion may 
contain odd, as well as even, powers of v. Thus, for example, v"^ should be replaced 
with V In phrases such as "expansion In powers of f ^" . 

The following paragraph should be Inserted after the paragraph that Includes Eq. (2.6): 

"The above estimates for the probabilities of \QQg) Fock states apply If the spin 
state of the QQ pair Is the same as In the dominant \QQ) Fock state. If the spin state Is 

1 



different, we must replace gA-'V in (2.6) with gBcr to obtain a nonzero matrix element. 
Using the velocity-scaling rules of Table I, we again obtain an estimate AE ~ Mv^ for 
the energy shift, implying that the probability for a \QQg) state containing a gluon 
with momentum on the order of Mv is Pqq„ ~ v^- However, in the derivation of the 
velocity-scaling rules in Ref. [14], it was assumed that dynamical gluons have momenta 
of order Mv. If the gluon has a much smaller momentum k, then the estimate M^f ^ 
for the operator gB in Table I should be replaced with k'^v'^. Using this to estimate 
the energy shift from a \QQg) Fock state containing a gluon with momentum of order 
Mf^, we obtain AE ~ Mv^ and Pqq„ ~ v^. Thus, gluons with very low momenta 
exhibit the suppression that is characteristic of the multipole expansion. We conclude 
that a \QQg) Fock state that can be reached from the dominant \QQ) Fock state by a 
spin-flip transition is dominated by dynamical gluons with momenta of order Mv and 
that the probability of such a Fock state is -Pgog ~ v^" 

The following paragraph should be added at the end of Sec. IID: 

"The above discussion applies to Fock states \QQg) in which the QQ pair has the 
same total spin quantum number S as in the dominant \QQ) state. The probabilities 
for Fock states \QQg) that can be reached from the dominant Fock state by a spin-flip 
transition also scale in a definite way with v. The probability for such a Fock state to 
contain a dynamical gluon with momentum of order Mv is of order f ^, just as in the 
case of a non-spin-fiip transition. However, in the case of a spin-flip transition, this 
momentum region dominates because, as we have seen, gluons with softer momenta, 
on the order of Mv ^, are suppressed by the multipole expansion. Thus, if the QQ pair 
in the dominant Fock state has angular-momentum quantum numbers '^^^^Lj, then 
the Fock state \QQg), with the QQ pair in a color-octet state with the same value of L 
but different total spin quantum number, has a probability of order v^. For example, 
if the dominant Fock state consists of a QQ pair in a ^5*1 state, then the Fock state 
\QQg) with the QQ pair in a color-octet ^5*0 state has a probability of order v^. If the 



dominant Fock state consists of a QQ pair in a ^Pi state, then the Fock state \QQg) 
with the QQ pair in a color-octet ^Pj state has probabihty of order v^." 

In the first paragraph of Sec. IIIA, the following two sentences should be inserted just 
before the last sentence of the paragraph: 

"The matrix element is suppressed by v^ relative to the velocity-scaling rules in Table I 
if On annihilates and creates QQ pairs in the same color-spin-orbital state as appears 
in one of the Fock states \QQg) that can be obtained from the dominant Fock state 
by a spin-flip transition. In such a Fock state, the QQ pair must be in a color-octet 
state with the same orbital-angular-momentum quantum number L as in the dominant 
\QQ) state, but with different total spin quantum number." 

After the flrst paragraph of Sec. IIIA, the following new paragraph should be inserted: 

"If perturbation theory remained accurate down to the scale Mv, then the spin- 
flip matrix elements would be suppressed by an additional power of v. The reason for 
this is that the contribution to a spin-flip matrix element that is suppressed by only 
v^ relative to the velocity-scaling rules is power ultraviolet divergent. Therefore, one 
could carry out a renormalization of the matrix element in which this contribution 
is subtracted. The corresponding contribution to the decay rate would then reside 
in the short-distance coefficient of the matrix element that is associated with the 
dominant Fock state. (Such a subtraction is carried out automatically if dimensional 
regularization is used to cut off the ultraviolet divergences in the matrix element.) 
Once the subtraction has been made, the leading contribution to the spin-ffip matrix 
element comes from the scale Mv'^. It is subject to the usual multipole suppression 
and scales as f ^ relative to the velocity-scaling rules. In practice, one usually makes 
such subtractions perturbatively. It is not clear, in the charmonium and bottomonium 
systems, that perturbation theory is sufficiently accurate at the scale Mv to remove 
the v^ contribution completely. Therefore, we assume in the error estimates below 



that the spin-flip matrix elements scale as v^ relative to the velocity-scaling rules." 

In the second paragraph of Sec. IllA, f^ should be replaced with v^ in the phrase 
"suppressed by f^ or more". In Eq. (3.1), the error estimate O(f^r) should be replaced 
with O(f^r). In the third paragraph of Sec. IIIA, f^ should be replaced with v^ in the 
phrase "are of order v'^T or higher" . 

In Eqs. (4.1a), (4.1b), (4.3a), and (4.3b), the error estimates should be 0{v^T). At 
the end of the paragraph containing Eq. (4.2), "relative order f^" should be replaced 
with "relative order v^" . 

In Eqs. (6.8a), (6.8b), (6.9a), and (6.9b), the error estimates should be 0{v^a). 

There is a typesetting error in Eqs. (3.19a) and (3.19b). The first factor on the 



right side should be J3Nc/27r, just as in Eqs. (3.19c) and (3.19d). In the subsequent 
sentence, "order f^" should be replaced with "relative order v"^" . 

In Eq. (5.4), the last color matrix should be T^ji. In Eq. (5.5), the coefficient of the 
second term on the right-hand side should be 4:/{N^ — 1), rather than 2/{N^ — 1). 

In Eqs. (A16) and (A25), the running coupling constant should be as{2M) rather than 
asiM). 
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